{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A finite difference model as a Python function\n", "*Prof. dr.ir.T.N.Olsthoorn*\n", "\n", "*Heemstede, Sept. 2016, May 2017*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generalize the finite difference model into a callable function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 3D finite difference computation in the previous chapter can be generalized by putting the generic parts in a callable Python function and turning the the user or case specific data into the arguments by which the function is called like so:\n", "\n", "Phi = fdm3( x, y, z, kx, ky, kz, FQ, IH, IBOUND )\n", "\n", "Doing so and adding some input error checking for convenience leads to the following code, that is save to disk for future use:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting fdm_a.py\n" ] } ], "source": [ "%%writefile fdm_a.py\n", "# write the function in this cell to disk as file fdm.py\n", "\n", "import numpy as np\n", "#import pdb # in case we need to debug this function\n", "\n", "def fdm3(x, y, z, kx, ky, kz, FQ, HI, IBOUND):\n", " '''Returns computed heads of steady state 3D finite difference grid.\n", " \n", " Steady state 3D Finite Difference Model that computes the heads a 3D ndarray.\n", " \n", " Parameters\n", " ----------\n", " `x` : ndarray, shape: Nx+1, [L]\n", " `x` coordinates of grid lines perpendicular to rows, len is Nx+1\n", " `y` : ndarray, shape: Ny+1, [L]\n", " `y` coordinates of grid lines along perpendicular to columns, len is Ny+1\n", " `z` : ndarray, shape: Nz+1, [L]\n", " `z` coordinates of layers tops and bottoms, len = Nz+1\n", " `kx`, `ky`, `kz` : ndarray, shape: (Ny, Nx, Nz) [L/T]\n", " hydraulic conductivities along the three axes, 3D arrays.\n", " `FQ` : ndarray, shape: (Ny, Nx, Nz), [L3/T]\n", " prescrived cell flows (injection positive, zero of no inflow/outflow)\n", " `IH` : ndarray, shape: (Ny, Nx, Nz), [L]\n", " initial heads. `IH` has the prescribed heads for the cells with prescribed head.\n", " `IBOUND` : ndarray of int, shape: (Ny, Nx, Nz), dim: [-]\n", " boundary array like in MODFLOW with values denoting\n", " * IBOUND>0 the head in the corresponding cells will be computed\n", " * IBOUND=0 cells are inactive, will be given value NaN\n", " * IBOUND<0 coresponding cells have prescribed head\n", " \n", " Returns\n", " ------- \n", " `Phi` : ndarray, shape: (Ny, Nx, Nz), [L]\n", " the 3D array with the final heads with `NaN` at inactive cells.\n", " \n", " TO 160905\n", " '''\n", "\n", " import numpy as np\n", " import scipy.sparse as sp\n", " from scipy.sparse.linalg import spsolve # to use its short name\n", "# pdb.set_trace()\n", " x = np.sort(np.array(x)) # enforce ascending\n", " y = np.sort(np.array(y))[::-1] # enforce descending\n", " z = np.sort(np.array(z))[::-1] # enforce descending\n", " \n", " # as well as the number of cells along the three axes\n", " SHP = Nz, Ny, Nx = len(z) - 1, len(y)-1, len(x)-1\n", "\n", " Nod = np.prod(SHP)\n", " \n", " if Nod == 0:\n", " raise AssertationError(\"Nx, Ny and Nz must be >= 1\")\n", "\n", " # assert correct shape of input arrays\n", " if kx.shape != SHP:\n", " raise AssertionError(\"shape of kx {0} differs from that of model {1}\".format(kx.shape,SHP))\n", " if ky.shape != SHP:\n", " raise AssertionError(\"shape of ky {0} differs from that of model {1}\".format(ky.shape,SHP))\n", " if kz.shape != SHP:\n", " raise AssertionError(\"shape of kz {0} differs from that of model {1}\".format(kz.shape,SHP))\n", " \n", " # from this we have the width of columns, rows and layers\n", " dx = np.abs(np.diff(x).reshape(1, 1, Nx)) # enforce positive\n", " dy = np.abs(np.diff(y).reshape(1, Ny, 1)) # enforce positive\n", " dz = np.abs(np.diff(z)).reshape(Nz, 1, 1) # enforce positive\n", " \n", " active = (IBOUND >0).reshape(Nod,) # boolean vector denoting the active cells\n", " inact = (IBOUND==0).reshape(Nod,) # boolean vector denoting inacive cells\n", " fxhd = (IBOUND <0).reshape(Nod,) # boolean vector denoting fixed-head cells\n", "\n", " # half cell flow resistances\n", " Rx = 0.5 * dx / (dy * dz) / kx\n", " Ry = 0.5 * dy / (dz * dx) / ky\n", " Rz = 0.5 * dz / (dx * dy) / kz\n", " \n", " # set flow resistance in inactive cells to infinite\n", " Rx = Rx.reshape(Nod,); Rx[inact] = np.Inf; Rx=Rx.reshape(SHP)\n", " Ry = Ry.reshape(Nod,); Ry[inact] = np.Inf; Ry=Ry.reshape(SHP)\n", " Rz = Rz.reshape(Nod,); Rz[inact] = np.Inf; Rz=Rz.reshape(SHP)\n", " \n", " # conductances between adjacent cells\n", " Cx = 1 / (Rx[:, :, :-1] + Rx[:, :, 1:])\n", " Cy = 1 / (Ry[:, :-1, :] + Ry[:, 1:, :])\n", " Cz = 1 / (Rz[:-1, :, :] + Rz[1:, :, :])\n", " \n", " NOD = np.arange(Nod).reshape(SHP)\n", " \n", " IE = NOD[:, :, 1:] # east neighbor cell numbers\n", " IW = NOD[:, :, :-1] # west neighbor cell numbers\n", " IN = NOD[:, :-1, :] # north neighbor cell numbers\n", " IS = NOD[:, 1:, :] # south neighbor cell numbers\n", " IT = NOD[:-1, :, :] # top neighbor cell numbers\n", " IB = NOD[1:, :, :] # bottom neighbor cell numbers\n", " \n", " R = lambda x : x.ravel() # generate anonymous function R(x) as shorthand for x.ravel()\n", "\n", " # notice the call csc_matrix( (data, (rowind, coind) ), (M,N)) tuple within tupple\n", " A = sp.csc_matrix(( -np.concatenate(( R(Cx), R(Cx), R(Cy), R(Cy), R(Cz), R(Cz)) ),\\\n", " (np.concatenate(( R(IE), R(IW), R(IN), R(IS), R(IB), R(IT)) ),\\\n", " np.concatenate(( R(IW), R(IE), R(IS), R(IN), R(IT), R(IB)) ),\\\n", " )),(Nod,Nod))\n", " \n", " # to use the vector of diagonal values in a call of sp.diags() we need to have it aa a \n", " # standard nondimensional numpy vector.\n", " # To get this:\n", " # - first turn the matrix obtained by A.sum(axis=1) into a np.array by np.array( .. )\n", " # - then take the whole column to loose the array orientation (to get a dimensionless numpy vector)\n", " adiag = np.array(-A.sum(axis=1))[:,0]\n", " \n", " Adiag = sp.diags(adiag) # diagonal matrix with a[i,i]\n", " \n", " RHS = FQ.reshape(Nod,1) - A[:,fxhd].dot(HI.reshape(Nod,1)[fxhd]) # Right-hand side vector\n", " \n", " Phi = HI.flatten() # allocate space to store heads\n", " \n", " Phi[active] = spsolve( (A + Adiag)[active][:,active] ,RHS[active] ) # solve heads at active locations\n", " \n", " Phi[inact] = np.NaN # put NaN at inactive locations\n", " \n", " return Phi.reshape(SHP) # reshape vector to 3D size of original model" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np # we always need this\n", "import fdm_a # import fdm_a module for use\n", "\n", "from importlib import reload # we need reload if we edited the file fdm_a.py\n", "reload(fdm_a) # if edited, must reload it\n", "\n", "fdm3 = fdm_a.fdm3 # for convenience create a local name to just write fdm3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This function should be saved in a .py file e.g. \"fdm3.py\" so that it can be used as a module that can be imported by other users or programs or scripts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Apply the model\n", "\n", "To apply the model we specify coofdinates, conductivities, prescribed flows and prescribed heads and the `IBOUND` array. Then we call the function with the proper arguments. After the function computed the heads, we'll show them by a contour plot.\n", "\n", "### Generate input to run the model with\n", "\n", "The lines below are case specific and could be placed in a script that can be run to set up the model afterw which the model is called from the script with the proper arguments" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# specify a rectangular grid\n", "x = np.arange(-1000., 1000., 25.)\n", "y = np.arange(-1000., 1000., 25.) # backward, i.e. first row grid line has highest y\n", "z = np.arange(-100., 0., 20.) # backward, i.e. from top to bottom\n", "\n", "SHP = (len(z)-1, len(y)-1, len(x)-1)\n", "\n", "k = 10.0 # m/d uniform conductivity\n", "kx = k * np.ones(SHP) # [L/T] 3D kx array\n", "ky = k * np.ones(SHP) # [L/T] 3D ky array with same values as kx \n", "kz = k * np.ones(SHP) # [L/T] 3D kz array with same values as kx\n", "\n", "FQ = np.zeros(SHP) # all flows zero. Note sz is the shape of the model grid\n", "FQ[2, 30, 25] = -1200 # [m3/d] extraction in this cell\n", "\n", "HI = np.zeros(SHP) # initial heads\n", "\n", "IBOUND = np.ones(SHP)\n", "IBOUND[:, -1, :] = -1 # last row of model heads are prescribed\n", "IBOUND[:, 40:45, 20:70]=0 # these cells are inactive" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Call the function with the correct arguments" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Phi = fdm3( x, y, z, kx, ky, kz, FQ, HI, IBOUND)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualization of the results: plot heads as contours" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import the required plotting module and setup the plot and run ``%matplotlib notebook`` to allow plots beeing shown inside the notebook." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib notebook\n", "import matplotlib.pyplot as plt # combines namespace of numpy and pyplot" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "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", " };\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 = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%matplotlib notebook\n", "import matplotlib.pyplot as plt # combines namespace of numpy and pyplot\n", "\n", "plt.figure()\n", "plt.xlabel('x [m]')\n", "plt.ylabel('head [m]')\n", "plt.title('Head between two fixed boundaries and constant precipitation {0} m/d'.format(N))\n", "plt.plot(xm, Phi.ravel(), 'bo-', label='numeric')\n", "\n", "# now add the analytical solution\n", "fi = N/(2*kD) * ((L/2)**2 - x**2)\n", "plt.plot(x, fi, 'xr', label='analytic')\n", "plt.legend()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows that the analytical and numerical solutions match. Notice, however, that we have made the outer two cells very narrow, so that the location of the `x[0]` and `xm[0]` and `x[-1]` and `xm[-1]` are almost indistinguishable, so that the position of the boundary are essentially the same for the numerical and the analytical model. (Just make the outer two cells as wide as the other cells to see the difference.\n", "\n", "The 3D model was used here as a 1D model along the `x`-axis. It should be evident that our 3D numerical model is very flexible in modeling 1D, 2D as well as 3D situations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 2, semi-confined flow (mazure case)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assume the case where we have a water body with constant head $h=h_w$ that at `x=0` is in direct contact with an semi aquifer that is charaterized with constant transmissivity `kD` convered by a semi-confining layer with constant resistance `c` above which we have a head that is maintained at a constant level $h=h_p$. This problem was solved by Marzure around 1930 when he studied groundwater flow from the IJssel Lake to a new adjacent polder with a 5 m lower maintained water level.\n", "\n", "The analytical solution for the head in the regional aquifer, that is, below the confining layer is\n", "\n", "$ \\phi - h_p = (\\phi_0 - h_p) \\,\\, \\exp (-x / \\lambda) $\n", "\n", "with $\\lambda = \\sqrt {kD c}$\n", "\n", "Let's work out this case to compare the analytical result with our model." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# import numpy as np, assumed numpy has been loaded above\n", "\n", "k = 25 # m/d conductivity of regional aquifer\n", "D = 50 # m thickness of regional aquifer\n", "d = 10 # m thickness of confining layer\n", "dtop = 0.001 # m dummy thickness of layer in which the head is maintained\n", "c = 500.0 # d vertical hydraulic resistance of confining layer\n", "hw = -0.4 # m fixed IJssel lake elevation\n", "hp = -5.0 # m fixed polder water elevation.\n", "\n", "kD = k*D\n", "lam = np.sqrt(kD * c)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "x = np.hstack((-0.001, np.linspace(0, 5000, 100)))\n", "y = np.array([-0.5, 0.5])\n", "grs = 0 # m, ground surface elevation\n", "z = np.array([grs-D-d-dtop, grs-d-dtop, grs-dtop, grs])" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "Nz, Ny, Nx = SHP = len(z)-1, len(y)-1, len(x)-1\n", "IBOUND = np.ones(SHP)\n", "IBOUND[0, :, :] = -1 # fixed head maintained in top layer (above confining layer)\n", "IBOUND[-1, :, 0] = -1 # fixed head in aquifer at x=0\n", "\n", "HI = np.zeros(SHP)\n", "HI[-1, :, 0] = hw # fixed level of IJssel Lake\n", "HI[ 0, :, :] = hp # fixed polder level\n", "\n", "FQ = np.zeros(SHP) # no fixed flows\n", "\n", "kAquif = k;\n", "kConf = d/c;\n", "kTop = 100; #immaterial\n", "\n", "K = np.array([kTop, kConf, kAquif]).reshape(Nz, 1, 1) * np.ones(SHP)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "Phi = fdm3(x, y, z, K, K, K, FQ, HI, IBOUND)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "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", " };\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 = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# %matplotlib notebook # already done above\n", "# import matplotlib.pyplot as plt # combines namespace of numpy and pyplot\n", "\n", "plt.figure()\n", "plt.xlabel('x [m]')\n", "plt.ylabel('head [m]')\n", "plt.title('Head between two fixed boundaries and constant precipitation {0} m/d'.format(N))\n", "\n", "xm = 0.5 * (x[:-1] + x[1:])\n", "plt.plot(xm, Phi[-1][0], 'bo-')\n", "\n", "# now add the analytical solution\n", "fi = hp + (hw - hp) * np.exp(-x/lam)\n", "plt.plot(x, fi, 'xr-')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As can be seen, we get same result as before, but now with only two layers, of which the top one represents the semi confined layer in which center the head is prescribed." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Cicular island with recharge\n", "\n", "This examples uses a flat model (one layer) and places a fixed head outside the boundary of the island." ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "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", " };\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 = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.linspace(-1000., +1000., 200)\n", "y = np.linspace(-1000., +1000., 200)\n", "z = np.array([0., -100.])\n", "\n", "xm = 0.5 * (x[:-1] + x[1:])\n", "ym = 0.5 * (y[:-1] + y[1:])\n", "zm = 0.5 * (z[:-1] + z[1:])\n", "\n", "Nx = len(xm)\n", "Ny = len(ym)\n", "Nz = len(zm)\n", "\n", "SHP = (Nz, Ny, Nx)\n", "\n", "ZM, YM, XM = np.meshgrid(zm, ym, xm, indexing='ij') # full 3D arrays of cell center coordinates\n", "\n", "DX = np.abs(np.diff(x).reshape((1, 1, Nx)) * np.ones(SHP)) # column width (3D aray)\n", "DY = np.abs(np.diff(y).reshape((1, Ny, 1)) * np.ones(SHP)) # row widths (3D array)\n", "\n", "x0 = 0.; y0 = 0. # center of the island\n", "\n", "RM = np.sqrt((XM - x0)**2 + (YM - y0)**2).reshape(SHP) # distance to center\n", "\n", "R = 750.0 # [m] radius of the island\n", "\n", "IBOUND = np.ones(SHP)\n", "IBOUND[RM>R] = -1\n", "\n", "k0 = 10 # [m/d]\n", "k = k0 * np.ones(SHP) # uniform conductivity\n", "\n", "kD = float(k0 * np.abs(np.diff(z)))\n", "\n", "rch = 0.01 # [m/d] recharge rate\n", "FQ = rch * DX * DY # [m3/d] cell inflows\n", "IH = np.zeros(SHP) # [m] initial heads\n", "\n", "Phi = fdm3(x, y, z, k, k, k, FQ, IH, IBOUND) # run model, return heads\n", "\n", "# plot the heads a contours\n", "plt.figure()\n", "plt.xlabel('x [m]'); plt.ylabel('y [m]')\n", "plt.title('Circular island with radius R={0} m and recharge N = {1} m/d'.format(R,N))\n", "plt.contourf(xm,ym,Phi[0],50)\n", "\n", "# plot the heads along the cross section and compare with analytical solution\n", "plt.figure()\n", "plt.xlabel('x [m]'); plt.ylabel('y [m]')\n", "plt.title('Heads through the center of the island')\n", "centerRow = int(np.floor(Ny/2))\n", "plt.plot(xm, Phi[0, centerRow, :], label='Numeric')\n", "Island = np.logical_and(xm >= -R, xm <= R)\n", "plt.plot(xm[Island], N/(4 * kD) * (R**2 - xm[Island] ** 2), label='Analytic')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Notice the difference between the analytic and numeric solution. This difference depends on the size of the number of cells that are used. The finer the model, the better the agreement. A small agreement remains because the boundary of the numeric model is not perfectly circular. We will a develop and use an axisymmetric model in the next chapter." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Circular polder\n", "This circular polder is also a circular island with fixed boundary at distance R from its center, but instead of a uniform recharge we have leakage through a semi-confining top layer that we will simulate using an extra layer.\n", "\n", "Within the island we have a fixed head above the confining top layer. As an extra example we can choose to extend the polder to not be surrounded by open water with a fixed head but by other land where the head is fixed above the confing layer but at a different value.\n", "\n", "Much of the set-up below is copied from the previous example." ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "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", " };\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 = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import scipy # we need bessel functions\n", "\n", "I0 = scipy.special.i0 # bessel function\n", "I1 = scipy.special.i1 # same\n", "K0 = scipy.special.k0 # same\n", "K1 = scipy.special.k1 # same\n", "\n", "z0= 0. # m, reference elevation, ground surface\n", "d = 10. # m, thickness of confining layer\n", "D = 100. # m, thickness of regional aquifer\n", "\n", "R = 750. # m, radius of island\n", "x0 = 0.; y0 = 0. # center of the island\n", "\n", "hp = -2.0 # m, maintained head r<=R\n", "hl = 0.50 # m, maintained head r>R\n", "\n", "c = 250. # d, vertical hydraulic resistance of confining layer\n", "k0 = d/c # m/d vertical conductivity of confining layer\n", "k1 = 10. # m/d, conductivity of regional layer\n", "kD = k1 * D # m2/d transmissivity of regional aquifer\n", "lam = np.sqrt(kD * c) # m, spreading or characteristic length of aquifer system\n", "\n", "# The grid, coordinates and size\n", "x = np.linspace(-2000., +2000., 200)\n", "y = np.linspace(-2000., +2000., 200)\n", "z = z0 - np.array([0, d, d+D])\n", "\n", "xm = 0.5 * (x[:-1] + x[1:])\n", "ym = 0.5 * (y[:-1] + y[1:])\n", "zm = 0.5 * (z[:-1] + z[1:])\n", "\n", "Nx = len(xm)\n", "Ny = len(ym)\n", "Nz = len(zm)\n", "\n", "SHP = (Nz, Ny, Nx)\n", "\n", "ZM, YM, XM = np.meshgrid(zm, ym, xm, indexing='ij') # full 3D arrays of cell center coordinates\n", "\n", "DX = np.abs(np.diff(x).reshape((1, 1, Nx)) * np.ones(SHP)) # column width (3D aray)\n", "DY = np.abs(np.diff(y).reshape((1, Ny, 1)) * np.ones(SHP)) # row widths (3D array)\n", "\n", "RM = np.sqrt((XM - x0)**2 + (YM - y0)**2).reshape(SHP) # distance to center\n", "in_polder = RM[0, :, :]<=R\n", "in_land = np.logical_not(in_polder)\n", "\n", "IBOUND = np.ones(SHP)\n", "IBOUND[0, :, :] = -1 # maintain head everywhere in confining layer\n", "\n", "k = np.ones(SHP) # uniform conductivity\n", "k[ 0, :, :] = k0/2 # m/d, top layer\n", "k[-1, :, :] = k1 # m/d, bottom layer\n", "\n", "rch = 0.01 # [m/d] recharge rate\n", "FQ = np.zeros(SHP)# [m3/d] cell inflows\n", "\n", "IH = np.zeros(SHP) # [m] initial heads\n", "H = np.ones((Ny, Nx)) # m, maintained head in top layer\n", "H[in_land] = hl # outside R\n", "H[in_polder]= hp # inside R\n", "\n", "IH[0, :, :] = H # use in initial heads\n", "\n", "Phi = fdm3(x, y, z, k, k, k, FQ, IH, IBOUND) # run model, return heads\n", "\n", "# plot the heads a contours\n", "plt.figure()\n", "plt.xlabel('x [m]'); plt.ylabel('y [m]')\n", "plt.title('Circular polder with surrounding land, semi confined regional aquifer')\n", "plt.contourf(xm,ym,Phi[-1],50) # contour bottom layer\n", "\n", "# plot the heads along the cross section and compare with analytical solution\n", "plt.figure()\n", "plt.xlabel('x [m]'); plt.ylabel('y [m]')\n", "plt.title('Heads through the center of the polder and land')\n", "centerRow = int(np.floor(Ny/2))\n", "\n", "plt.plot(xm, Phi[-1, centerRow, :], label='Numeric')\n", "Island = np.logical_and(xm>=-R, xm<=R)\n", "\n", "# Analytical solution\n", "RL = R/lam # shorthand\n", "\n", "# head at x==R\n", "xPold = xm[np.abs(xm) <= R] - x0 # inside polder\n", "xLand = xm[np.abs(xm) > R] - y0 # outside polder\n", "\n", "# head a x==R, analytical\n", "phiR = (hp * I1(RL) / I0(RL) + hl* K1(RL) / K0(RL)) / (I1(RL) / I0(RL) + K1(RL) / K0(RL))\n", "\n", "phiPold = hp + (phiR - hp) * I0(np.abs(xPold) / lam) / I0(RL) # h in polder, r<=R\n", "phiLand = hl + (phiR - hl) * K0(np.abs(xLand) / lam) / K0(RL) # h outside polder, r> R\n", "\n", "plt.plot(xPold, phiPold, 'r-', label='Poldm analytic')\n", "plt.plot(xLand, phiLand, 'g-', label='Land analytic')\n", "\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The analytical solution matches the numerical one except near the edges of the model, where the model is closed and the analytical solution is not. The numerical solution, therefore, becomes horizontal near the boundary, indicating no flow across it. One also notices the distortion of the circular color field in the first figure due to the boundaries of the square shape of our model. Notice that the only boundary conditions in this model are prescribed heads in the top layer, which are different inside and outside the radius `R`.\n", "\n", "One nicely see both analytical solutions, one for the land outside the polder (green) and one for the land inside teh polder (red).\n", "\n", "There is a small difference between the numerical and analytical solution in the center. Again, this may be due to the fact that the circular polder can't be perfectly circular due to the rectangular cells of the model. But the larger the number of cells, the smaller the difference will be.\n", "\n", "It is important to note that we used half the vertical conductivity of the top layer. This is because in the model, water flows vertically from the center of the toplayer, where the head is fixed to the bottom, that is, it passes only half the cell. To keep the resistance of this lower half of the cell equal to the given vertical hydraulic resisttance, we have to use $ k0 = (d/2) / c $ instead of $ k0 = d/c $" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## More efficient coordinates\n", "One sees that to set up a model, there area quite some lines that deal with the grid coordinates. This can be done a lot more efficient with a Grid class, the instances of which are invoked with the grid coordinates, after which a large number of variables can be requested from it, variables that are computed upon request." ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [conda root]", "language": "python", "name": "conda-root-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }