{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Computing flows with the finite difference model\n", "\n", "*Prof. dr.ir.T.N.Olsthoorn*\n", "\n", "*Heemstede, Sept. 2016*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding flows to the output of fdm3\n", "\n", "Until now, the 3D finite difference model was used to compute only the heads in the cell centers. It should, however been clear that the model has all the information on board to also compute the flows. In fact, each line in the system matrix is a cell's water balance and computes flows over all faces of the cell in question including the total net inflow of each cell. The equation\n", "\n", "$ Q_x = -C_x \\times \\,\\, diff (Phi, axis=1) $\n", "\n", "is the flow across all cell faces perpendicular to the `x-axis`. Likewise\n", "\n", "$ Q_y = +C_y \\times \\, \\, diff (Phi, axis=0) $\n", "\n", "is the flow across all cell faces perpendicular to the `y-axis` in direction of this axis, and\n", "\n", "$ Q_z = +C_z \\times \\, \\, diff(Phi, axis=2) $\n", "\n", "is the upward vertical flow\n", "\n", "In fact, we can readily compute these flow across cell faces within our model as is shown below.\n", "\n", "We can also compute the total net inflow of all individual cells from the matrix multiplication\n", "\n", "$ Q = A \\times \\Phi $\n", "\n", "where `A` is the system matrix with the conductances and $\\Phi$ is the complete vector of heads, including the fixed heads, as all computed and prescribed heads are known after the model has been solved.\n", "\n", "## Flow output of the model\n", "\n", "Because we have all the information to compute all flows at our disposition inside the function when computing the heads, we can, and perhaps should, at the same time compute these important flow arrays:\n", "\n", "* `Q` : ndarray shape (`Ny`, `Nx` , `Nz`) [L3/T] # total net inflow of cells\n", "* `Qx` : ndarray shape (`Ny`, `Nx`-1, `Nz`) [L3/T] # flow across cell faces in `x`-direction\n", "* `Qy` : ndarray shape (`Ny`-1, `Nx`, `Nz`) [L3/T] # flow across cell faces in `y`-direction\n", "* `Qz` : ndarray shape (`Ny`, `Nx`, `Nz`-1) [L3/T] # flow across cell faces in `z`-direction\n", "\n", "We gather the computed arrays of `Phi`, `Q`, `Qx`, `Qy` and `Qz` in an single `named tuple` output named `Out`, from which the individual arrays may be addressed by their name like so:\n", "\n", "`Out.Phi`, `Out.Q`, `Out.Qx`, `Out.Qy` and `Out.Qz`\n", "\n", "### Flows in the cell centers\n", "\n", "The flows at across the cell iterfaces and the net inflow for each cell is uniquely defined by the equations given above. But this is not true for the velocity or specific discharge, because the cell face area jumps at the cell itnerface. So we mat have a different specific discharge perpendicular to the cell face just to the right and to the left of it.\n", "\n", "Therefore, such specific discharge (the Darcy velocity) are generally computed for the cell centers by averaging the total flow at two opposite cell faces and deviding by the cross section perpendicular to the cell face. This vector is unique, yet not completely defined for the case when we also have an extraction in the cell. But this does not geneally or necessarly disturb the picture of velocity vectors much. So we will use this method further down.\n", "\n", "We have written a function `quivdata` to extract the necessary coordinates and flows or velocities for displaying them as arrows with the matplotlib function `quiver`, see example further down.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Some more additions been made to the model:\n", "\n", "- Applying the function `unique` to `x`, `y` and `z` input vectors to prevent duplicates. The function unique is defined in the module `fdm`. It allows entering coordinates in arbitrary order as they will be sorted and double values will be eliminated according to some user definable tolerance.\n", "\n", "- Keeping track of the axis directions. We demand that `x` must be ascending, while `y` and `z` must be descending. This way columns (x) run from left to right when printing, the top row is the one with the largest `y`-coordinate and the top layer is the one with the highest `z`-coordinate. Enforcing these directions make interpreting scrolled or printed output intuitive and ensures the correct sign of flows. That is, a postive `Qx`, `Qy` and `Qz` always points into the direction of higher coordinates.\n", "\n", "The adapted model is shown below.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## fdm3 model with computation of flows\n", "\n", "The model `fdm3` is written to the module `fdm.py` together with the functions `quivdata` and `unique`" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting fdm_b.py\n" ] } ], "source": [ "%%writefile fdm_b.py\n", "\n", "import numpy as np\n", "import pdb\n", "import scipy.sparse as sp\n", "from scipy.sparse.linalg import spsolve # to use its short name\n", "from collections import namedtuple\n", "\n", "class InputError(Exception):\n", " pass\n", "\n", "def quivdata(Out, x, y, z=None, iz=0):\n", " \"\"\"Returns coordinates and velocity components to show velocity with quiver\n", " \n", " Compute arrays to display velocity vectors using matplotlib's quiver.\n", " The quiver is always drawn in the xy-plane for a specific layer. Default iz=0\n", " \n", " Parameters\n", " ----------\n", " `Out` : namedtuple holding arrays `Qx`, `Qy`, `Qz` as defined in `fdm3`\n", " `Qx` : ndarray, shape: (Ny, Nx-1, Nz), [L3/T]\n", " Interfacial flows in finite difference model in x-direction from `fdm3'\n", " `Qy` : ndarray, shape: (Ny-1, Nx, Nz), [L3/T]\n", " Interfacial flows in finite difference model in y-direction from `fdm3`\n", " `Qz` : ndarray, shape: (Ny, Nx, Nz-1), [L3/T]\n", " Interfacial flows in finite difference model in z-direction from `fdm3` \n", " `x` : ndarray, [m]\n", " Grid line coordinates of columns\n", " 'y' : ndarray, [m]\n", " Grid line coordinates of rows\n", " `z` : ndaray [L] | int [-]\n", " If z == None, then iz must be given (default = 0)\n", " If z is an ndarray vector of floats\n", " z will be interpreted as the elvations of uniform layers.\n", " iz will be ignored\n", " If z is a full 3D ndarray of floats\n", " z will be interpreted as the elevations of the tops and bottoms of all cells.\n", " iz will be ignored\n", " `iz` : int [-]\n", " iz is ignored if z ~= None\n", " iz is the number of the layer for which the data are requested,\n", " and all output arrays will be 2D for that layer.\n", "\n", " Returns\n", " -------\n", " `Xm` : ndarray, shape: (Nz, Ny, Nx), [L]\n", " x-coordinates of cell centers\n", " `Ym` : ndarray, shape: (Nz, Ny, Nx), [L]\n", " y-coodinates of cell centers\n", " `ZM` : ndarray, shape: (Nz, Ny, Nx), [L]\n", " `z`-coordinates at cell centers\n", " `U` : ndarray, shape: (Nz, Ny, Nx), [L3/d]\n", " Flow in `x`-direction at cell centers\n", " `V` : ndarray, shape: (Nz, Ny, Nx), [L3/T]\n", " Flow in `y`-direction at cell centers\n", " `W` : ndarray, shape: (Nz, Ny, Nx), [L3/T]\n", " Flow in `z`-direction at cell centers.\n", " \n", " \"\"\"\n", " Ny = len(y)-1\n", " Nx = len(x)-1\n", " \n", " xm = 0.5 * (x[:-1] + x[1:])\n", " ym = 0.5 * (y[:-1] + y[1:])\n", " \n", " X, Y = np.meshgrid(xm, ym) # coordinates of cell centers\n", " \n", " # Flows at cell centers\n", " U = np.concatenate((Out.Qx[iz, :, 0].reshape((1, Ny, 1)), \\\n", " 0.5 * (Out.Qx[iz, :, :-1].reshape((1, Ny, Nx-2)) +\\\n", " Out.Qx[iz, :, 1: ].reshape((1, Ny, Nx-2))), \\\n", " Out.Qx[iz, :, -1].reshape((1, Ny, 1))), axis=2).reshape((Ny,Nx))\n", " V = np.concatenate((Out.Qy[iz, 0, :].reshape((1, 1, Nx)), \\\n", " 0.5 * (Out.Qy[iz, :-1, :].reshape((1, Ny-2, Nx)) +\\\n", " Out.Qy[iz, 1:, :].reshape((1, Ny-2, Nx))), \\\n", " Out.Qy[iz, -1, :].reshape((1, 1, Nx))), axis=1).reshape((Ny,Nx))\n", " return X, Y, U, V\n", "\n", "\n", "def unique(x, tol=0.0001):\n", " \"\"\"return sorted unique values of x, keeping ascending or descending direction\"\"\"\n", " if x[0]>x[-1]: # vector is reversed\n", " x = np.sort(x)[::-1] # sort and reverse\n", " return x[np.hstack((np.diff(x) < -tol, True))]\n", " else:\n", " x = np.sort(x)\n", " return x[np.hstack((np.diff(x) > +tol, True))]\n", "\n", " \n", "def fdm3(x, y, z, kx, ky, kz, FQ, HI, IBOUND):\n", " '''Steady state 3D Finite Difference Model returning computed heads and flows.\n", " \n", " Heads and flows are returned as 3D arrays as specified under output parmeters.\n", " \n", " Parameters\n", " ----------\n", " `x` : ndarray,[L]\n", " `x` coordinates of grid lines perpendicular to rows, len is Nx+1\n", " `y` : ndarray, [L]\n", " `y` coordinates of grid lines along perpendicular to columns, len is Ny+1\n", " `z` : ndarray, [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, shape: (Ny, Nx, Nz) of int\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", " outputs\n", " ------- \n", " `Out` : namedtuple containing heads and flows:\n", " `Out.Phi` : ndarray, shape: (Ny, Nx, Nz), [L3/T] \n", " computed heads. Inactive cells will have NaNs\n", " `Out.Q` : ndarray, shape: (Ny, Nx, Nz), [L3/T]\n", " net inflow in all cells, inactive cells have 0\n", " `Out.Qx : ndarray, shape: (Ny, Nx-1, Nz), [L3/T] \n", " intercell flows in x-direction (parallel to the rows)\n", " `Out.Qy` : ndarray, shape: (Ny-1, Nx, Nz), [L3/T] \n", " intercell flows in y-direction (parallel to the columns)\n", " `Out.Qz` : ndarray, shape: (Ny, Nx, Nz-1), [L3/T] \n", " intercell flows in z-direction (vertially upward postitive)\n", " the 3D array with the final heads with `NaN` at inactive cells.\n", " \n", " TO 160905\n", " '''\n", "\n", " # define the named tuple to hold all the output of the model fdm3\n", " Out = namedtuple('Out',['Phi', 'Q', 'Qx', 'Qy', 'Qz'])\n", " Out.__doc__ = \"\"\"fdm3 output, , containing fields Phi, Qx, Qy and Qz\\n \\\n", " Use Out.Phi, Out.Q, Out.Qx, Out.Qy and Out.Qz\"\"\" \n", "\n", " x = unique(x)\n", " y = unique(y)[::-1] # unique and descending\n", " z = unique(z)[::-1] # unique and 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 AssetError(\n", " \"Grid shape is (Ny, Nx, Nz) = {0}. Number of cells in all 3 direction must all be > 0\".format(SHP))\n", " \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)\n", " dy = np.abs(np.diff(y)).reshape(1, Ny, 1)\n", " dz = np.abs(np.diff(z)).reshape(Nz, 1, 1)\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", " # also notice that Cij = negative but that Cii will be postive, namely -sum(Cij)\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", " Out.Phi = HI.flatten() # allocate space to store heads\n", " \n", " Out.Phi[active] = spsolve( (A+Adiag)[active][:,active] ,RHS[active] ) # solve heads at active locations\n", " \n", " # net cell inflow\n", " Out.Q = (A+Adiag).dot(Out.Phi).reshape(SHP)\n", "\n", " # set inactive cells to NaN\n", " Out.Phi[inact] = np.NaN # put NaN at inactive locations\n", " \n", " # reshape Phi to shape of grid\n", " Out.Phi = Out.Phi.reshape(SHP)\n", " \n", " #Flows across cell faces\n", " Out.Qx = -np.diff( Out.Phi, axis=2) * Cx\n", " Out.Qy = +np.diff( Out.Phi, axis=1) * Cy\n", " Out.Qz = +np.diff( Out.Phi, axis=0) * Cz\n", " \n", " return Out # all outputs in a named tuple for easy access" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These functions have now been saved in the file \"fdm_b.py\". This file is a module, which can be imported, after which the functions in it can be used in programs and scripts." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import fdm_b\n", "from importlib import reload\n", "reload(fdm_b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Application of the model\n", "\n", "### Generate input to run the model with (same example)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "reload(fdm_b) # make sure we reload because we edit the file regularly\n", "\n", "# specify a rectangular grid\n", "x = np.arange(-1000., 1000., 25.)\n", "y = np.arange( 1000., -1000., -25.)\n", "z = np.array([20, 0 ,-10, -100.])\n", "\n", "Nz, Ny, Nx = 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)\n", "ky = k * np.ones(SHP)\n", "kz = k * np.ones(SHP)\n", "\n", "IBOUND = np.ones(SHP)\n", "IBOUND[:, -1, :] = -1 # last row of model heads are prescribed\n", "IBOUND[:, 40:45, 20:70]=0 # inactive\n", "\n", "FQ = np.zeros(SHP) # all flows zero. Note SHP is the shape of the model grid\n", "FQ[1, 30, 25] = -1200 # extraction in this cell\n", "\n", "HI = np.zeros(SHP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Call the function with the correct arguments" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Out = fdm_b.fdm3( x, y, z, kx, ky, kz, FQ, HI, IBOUND)\n", "\n", "# when not using some of the parameters\n", "#Phi, _, _, _ = fdm.fdm3( x, y, z, kx, ky, kz, FQ, HI, IBOUND) # ignore them using the _\n", "#Phi, _Qx, _Qy, _Qz = fdm.fdm3( x, y, z, kx, ky, kz, FQ, HI, IBOUND) # make them private\n" ] }, { "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." ] }, { "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": [ { "name": "stdout", "output_type": "stream", "text": [ "Out.Phi.shape = (3, 79, 79)\n", "Out.Q.shape = (3, 79, 79)\n", "Out.Qx.shape = (3, 79, 78)\n", "Out.Qy.shape = (3, 78, 79)\n", "Out.Qz.shape = (2, 79, 79)\n" ] } ], "source": [ "print('Out.Phi.shape = {0}'.format(Out.Phi.shape))\n", "print('Out.Q.shape = {0}'.format(Out.Q.shape))\n", "print('Out.Qx.shape = {0}'.format(Out.Qx.shape))\n", "print('Out.Qy.shape = {0}'.format(Out.Qy.shape))\n", "print('Out.Qz.shape = {0}'.format(Out.Qz.shape))" ] }, { "cell_type": "code", "execution_count": 7, "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": { "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": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plt.figure(); plt.title('Q of cells in top layer [m3/d]')\n", "plt.imshow(Out.Q[0, :, :], interpolation='None')\n", "plt.colorbar()\n", "\n", "plt.figure(); plt.title('Q of cells in second layer [m3/d]')\n", "plt.imshow(Out.Q[1, :, :], interpolation='None')\n", "plt.colorbar()\n", "\n", "plt.figure(); plt.title('Q of cells in bottom layer with well [m3/d]')\n", "plt.imshow(Out.Q[1, :, :], interpolation='None')\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice:\n", "\n", "* that the coordinates are the cell numbers, the model is (79, 79, 3).\n", "* the colorbars are different for each figure due to aut-selection.\n", "* the flows in the 3rd layer with the well varies indeed between -1200 ( well extraction) and 0\n", "* Further only the lower boundary is colored, i.e.non-zero, which is due to the fact that the heads are prescribed along that boundary and, therefore, the flows are computed. Positive values mean water is flowing from the outside world into the model.\n", "* The zone with inactive cells are not visible, because the flow in those cells are zero as is the case for other cells with no prescribed flows or heads." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Conclusion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have extended the 3D steady state finite difference model to not only compute the heads but also the flows between the cells and the net inflow of the cells. The latter are for the computed heads as well as for the prescribed heads. Inactive cells will have a head `NaN` (Not a Number) and will have flow zero.\n", "\n", "We also now have the convenience function to extract the flows at the center of all cells for visualization using the `matplotlib` function `quiver()`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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 }