{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Transient flow\n", "*Prof. dr.ir. T.N.Olsthoorn*\n", "\n", "*Heemstede, Oct. 2016, 24 May 2017*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Transient flow, in fact, requires a relatively straightforward extension of our steady-state model fdm3. The water balance equation for an arbitrary single cell in our model becomes:\n", "\n", "$$ \\intop_t^{t+\\Delta t}\\left( \\sum Q_{in} + Q_{ext} \\right) dt = W_{t+\\Delta t} - W_t $$ \n", "\n", "The the left side of the equation describes the total net inflow for this cell, divided in a sum that originates from the surrounding cells and the $ Q_{ext} $ denotes the total net inflow from the outside world (Injection, Recharge, Leakage etc.). It's clear that extractions are counted as negative injections. This total net inflow shoul be integrated over the considered time spand, i.e. $ t \\rightarrow t + \\Delta t $, during which the flows may vary in arbitrary ways. To the right we have the volume of water in the cell at the end of this period minus that at the start of this period. There is no approximation in this equation.\n", "\n", "When it comes to our model, we will write the right hand side in terms of volume, storage coefficient $ S_S $, and head $ h $. The integral to the left is dropped by writing for the flow their average values during the considered period:\n", "\n", "$$ \\left( \\sum \\overline{Q}_{in} + \\overline{Q}_{ext} \\right) \\Delta t \n", "= S_S V \\left( h_{t+\\Delta t} - h_t \\right)\n", "$$ \n", "\n", "where\n", "\n", "$$ V = \\Delta x \\Delta y \\Delta z $$\n", "\n", "the volume of the cell.\n", "\n", "As we saw in an earlier chapter where we derived the finite differenc method, the flow from each neighbor indexed $ j $ into the considered cell with index $ i $ is formulated as\n", "\n", "$$ Q_{ji} = C_{ij} (h_j - h_i) $$ \n", "\n", "Thus, $ Q $ varies during the time step according to how the heads vary. The question then is, which $ h $ is the true average during the time step, such that, then it is used, we get the average flows during the time step and hence we can integrate by multiplying with $ \\Delta t $.\n", "\n", "The answer is, we don't know. We only know for sure, that there exists some value of $ 0< \\epsilon < 1$ such that\n", "\n", "$$ Q_{t + \\epsilon \\Delta t } = \\overline{Q} $$\n", "\n", "But then, if we solve the system of equation, we wil end up with the head for this time, $ h_{t + \\epsilon \\Delta t} $, which definitely differs from the head at the end of the time step, which we need to compute the storage during the time step as expressed in the equation above.\n", "\n", "To solve this dilemma, we have to make an assumption, which is, that we assume that the time steps of our simulation are small enough to safely approximate the change of head during them as being linear. With this approximation, we have\n", "\n", "$$ h_{t+\\epsilon \\Delta t} = h_t + \\epsilon (h_{t + \\Delta t} - h_t) $$\n", "\n", "or\n", "\n", "$$ h_{t + \\Delta t} = h_t + \\frac { h_{t + \\epsilon \\Delta t} - h_t } \\epsilon $$\n", "\n", "This implies that we can now replace the unknown head at the end of the time step by the one we actually compute during the time step, by writing\n", "\n", "$$ \\left( \\sum \\overline{Q}_{in} + \\overline{Q}_{ext} \\right)_{t=t+\\epsilon \\Delta t} \n", "= S_S \\frac V {\\epsilon \\Delta t} \\left( h_{t+ \\epsilon \\Delta t} - h_t \\right)\n", "$$ \n", "\n", "If we just realize that all flows and heads will be compute that $ t + \\epsilon \\Delta t $, we can omit this index and write and reorder keeping the unknown $ Q_{in} $ to the left and putting all known $ Q_{ext} $ to the right, positive when entering the cell, we obtain:\n", "\n", "$$ - \\sum Q_{in} = Q_{ext} - S_S \\frac V {\\epsilon \\Delta t} \\left( h - h_t \\right)\n", "$$ \n", "\n", "We see that $ h_t $ is known, because it is the head in the considered cell at the beginning of the considered time step, i.e. at the end of the last time step. Therefore, we can rearrange to have all known values to the right of the equal sign and all unknowns at the left. This gives\n", "\n", "$$ - \\sum Q_{in} + S_S \\frac V {\\epsilon \\Delta t} h = Q_{ext} + S_S \\frac V {\\epsilon \\Delta t} h_t\n", "$$\n", "\n", "\n", "Just remember that the first term, $ -\\sum Q_{in} $ can be written out as a vector product, this becomes\n", "\n", "$$\n", "\\left[\n", "\\begin{array}{cccc}\n", "-C_E, & \\ldots & -C_B, & + \\left( C_E + \\ldots + C_B \\right)\n", "\\end{array}\n", "\\right]\n", "\\left[\n", "\\begin{array}{c}\n", "h_E \\\\ h_W \\\\ h_N \\\\ h_S \\\\ h_T \\\\ h_B \\\\ h\n", "\\end{array}\n", "\\right]\n", "+ S_S \\frac V {\\epsilon \\Delta t} h = Q_{ext} + S_S \\frac V {\\epsilon \\Delta t} h_t\n", "$$\n", "\n", "To bring the right-most term of the left-hand side of this equation under the vector product, we only have to put the the factor in front of $ h $ to the sum of coefficients like so:\n", "\n", "$$\n", "\\left[\n", "\\begin{array}{cccc}\n", "-C_E, & \\ldots & -C_B, & + \\left( C_E + \\ldots + C_B + S_S \\frac V {\\epsilon \\Delta t} \\right)\n", "\\end{array}\n", "\\right]\n", "\\left[\n", "\\begin{array}{c}\n", "h_E \\\\ h_W \\\\ h_N \\\\ h_S \\\\ h_T \\\\ h_B \\\\ h\n", "\\end{array}\n", "\\right]\n", " = Q_{ext} + S_S \\frac V {\\epsilon \\Delta t} h_t\n", "$$\n", "\n", "\n", "We see that we have now a water balance equation for the transient cell that has exactly the same structure as the one we developed for the steady-state case, with unknowns at the left and the known flows on the right. Notice that the second term to the right has the same dimension as $Q_{ext} $, namely [L3/T]. The include the storage due to the unknown head, it suffices to put the coefficient $ S_S V / (\\epsilon \\Delta t) $ in the coefficient of the requested node to the left. It doesn't change any of the other coefficients.\n", "\n", "If we compare this with the first chapter, in which it was shown how exchange between model and a fixed head in the outside world through resistance was implemented, we see that it's the same with implementation of transient flow. The conducntace is added to the matrix coefficient on the left, while the known part C_{ext}h_{ext} is kept on the right side, and it too has dimension [L3/T].\n", "\n", "Consedering the entire model, we have such an equation for each and every cell in it. These equations are combined into a system of equations\n", "\n", "$$ \\mathbf{A} \\times \\mathbf{h}_{t+\\epsilon \\Delta t} = \\overline{\\mathbf{Q}}_{ext} +\n", "\\mathbf{S}_S \\frac {\\mathbf{V}} {\\epsilon \\Delta t} \\mathbf{h}_t $$\n", "\n", "Where \\matbh{A} is the system matrix, in which the coefficient vector $ \\mathbf{S}_S \\mathbf{V}/(\\epsilon \\Delta t) $ was added to its diagonal.\n", "\n", "It is clear that these coefficents change with the lenght of the time step and, therefore the diagonal and the right-hand side of the equation have to be adapted with each new time step, if its length changes.\n", "\n", "The dealing with fixed heads does not change and was explained earlier in the chaper \"finite differnce modeling\".\n", "\n", "The answer that we obtain after having solved the system equation, is the head in all cells at time $ t + \\epsilon \\Delta t $. Therefore we a small extra step to obtain the head $ h_{t + \\Delta t} $ at the end of the time step\n", "\n", "$$ \\mathbf{h}_{t+\\Delta t} = \\mathbf{h}_t + \\frac { \\mathbf{h}_{t + \\epsilon \\Delta t} - \\mathbf{h}_t} \\epsilon $$\n", "\n", "\n", "The only point that we have still circumvented is the choice of $ 0\\epsilon<1 $, the so called implicitness. It is the point in time expressed as the fraction of the current time-step length at which the flows and heads represent the average values during the entire time step. In fact, we don't know. If the head change is linear, then $ \\epsilon=0.5 $ would be exact. But when the head exponentially appoaches a new equilibrium a value $ \\epsilon > 0.5 $ is better. Not only because heads will always thrive to a new equilibrium with time, but also because values of $ \\epsilon<0.5 $ yield unstable solutions, we always choose a value larger than 0.5. The most drastic choice is $\\epsilon = 1$, a choice which is said to make the model fully implicit. It conceptually assumes that the heads at the end of the current timestep well represent their average values during the time step. It may not be the most accurate value, especially with larger time steps, but it makes the model rock stable. In fact this is the choice that the world's most used finite difference model, MODFLOW, implicitly makes.\n", "\n", "Is it a good or a bad choice? Well it's not bad as errors due to this time discretization well damp out during subsequent time steps. Nevertheless, we will keep the value of $ \\epsilon $ explicit, so as to allow investigating the sensitivity of the outcomes of the model for it." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " ## Procedure\n", " \n", "We compute the required coefficients for each time step as follows:\n", "\n", "Compute\n", "\n", "for non-axially symmetric models:\n", "\n", "$$ \\mathbf{C}_S = \\frac {\\mathbf{S}_S \\mathbf{\\Delta x \\Delta y \\Delta z}} {\\epsilon \\Delta t} $$\n", "\n", "for axially symmetric models:\n", "\n", "$$ \\mathbf{C}_S = \\frac {\\mathbf{S}_S \\pi \\left( \\mathbf{x}[1:]^2 - \\mathbf{x}[:-1]^2 \\right) \\Delta z } {\\epsilon \\Delta t} $$\n", "\n", "Add this vector $ \\mathbf{C}_S $ explicitly to the diagonal of the system matrix at each time step.\n", "\n", "Compute\n", "\n", "$$ \\mathbf{C}_S \\mathbf{h}_t $$\n", "\n", "at the RHS of the system equation for each time step.\n", "\n", "We solve for $ \\mathbf{h}_{t + \\epsilon \\Delta t} $\n", "\n", "$$\n", "\\mathbf{h}_{t+\\epsilon \\Delta t} = \\mathbf{A}^{-1} \\times \\left( \\mathbf{Q}_{ext} +\n", "\\mathbf{C}_S \\mathbf{h}_t \\right)\n", "$$\n", "\n", "Finally compute\n", "\n", "$$ \\mathbf{h}_{t+\\Delta t} = \\mathbf{h}_t + \\frac {\\mathbf{h}_{t + \\epsilon \\Delta t} - \\mathbf{h}_t } \\epsilon $$\n", "\n", "And just before computing the new time step, update the know $\\mathbf{h}_t$ with the value just computed\n", "\n", "$$ \\mathbf{h}_t = \\mathbf{h}_{t+\\Delta_t} $$" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Implementation\n", "\n", "We have to adapt our fdm3 model at a few places. First is its signature. We need two extra inputs, namely $t$ and $ \\mathbf{S}_S and a possibility to adapt the implicitness $ \\epsilon $. Calling this model `fdm3t` we have\n", "\n", "Out = fdm3t(gr, t, kx, ky, kz, Ss, FQ, FH, IBOUND, epsilon=0.67)\n", "\n", "Where the output Out contains the computed arrays, see doc string of implementation below.\n", "\n", "The heads are at the end of the time steps not including the intial heads.\n", "The flows [L3/T] are average flows during each time step.\n", "\n", "All these array are, therefore, 4 dimensional, 3 spacial dimensions and and the fourth being time. The shape of the heads, $ Q $ and $ Qs $ arrays are, therore, $ (Ny, Nx, Nz, Nt) $\n", "\n", "The new $ Qs $ is the flow [m3/T] during the time step that enters the cell from storage (because all flows are positive when they enter a cell). Thus $Qs$ is positive when the head declines, yielding water from storage to the cell that then flows towards surrounding cells (or to the external world). \n", "\n", "We only use the specific storage of each cell. This can be easily computed from the total storage and even from the specific yield when required.\n", "\n", "What we have not (yet) implemented are non-lineairities like change of transmissivity of the model and, therefore, its cells under transient unconfined conditions. To do this is not difficult. It requires updating the transmissivities of the cells after a given number of so-called inner iterations. Each such adaptation is called an outer iterations. Once the head and or flow changes after an out iterations have become negligible, the model is said to have converged and the outcomes are used. Notice, however, that convergence is not always guaranteed for non-linear systems that are solved in terms of linear equations that are updated like it is done here as well as in MODFLOW. Special solvers can do a better job by adding non-linear Newton-Raphson schemes to the solver. MODFLOW.NWT and the new MODFLOW.USG have such a solver, which may be necessary for strongly non-linear models. These issues are beyond this course." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation; the adepted module to include axial symmetry" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting fdm_t.py\n" ] } ], "source": [ "%%writefile fdm_t.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 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 fdm3t(gr, t, kx, ky, kz, Ss, FQ, HI, IBOUND, epsilon=0.67):\n", " '''Transient 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", " 'gr' : `grid_object`, generated by gr = Grid(x, y, z, ..)\n", " if `gr.axial`==True, then the model is run in axially symmetric model\n", " t : ndarray, shape: [Nt+1]\n", " times at which the heads and flows are desired including the start time,\n", " which is usually zero, but can have any value.\n", " `kx`, `ky`, `kz` : ndarray, shape: (Nz, Ny, Nx), [L/T]\n", " hydraulic conductivities along the three axes, 3D arrays.\n", " `Ss` : ndarray, shape: (Nz, Ny, Nx), [L-1]\n", " specific elastic storage\n", " `FQ` : ndarray, shape: (Nz, Ny, Nx), [L3/T]\n", " prescrived cell flows (injection positive, zero of no inflow/outflow)\n", " `IH` : ndarray, shape: (Nz, Ny, Nx), [L]\n", " initial heads. `IH` has the prescribed heads for the cells with prescribed head.\n", " `IBOUND` : ndarray, shape: (Nz, Ny, Nx) 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", " `epsilon` : float, dimension [-]\n", " degree of implicitness, choose value between 0.5 and 1.0\n", " \n", " outputs\n", " ------- \n", " `Out` : namedtuple containing heads and flows:\n", " `Out.Phi` : ndarray, shape: (Nt+1, Nz, Ny, Nx), [L3/T] \n", " computed heads. Inactive cells will have NaNs\n", " To get heads at time t[i], use Out.Phi[i]\n", " Out.Phi[0] = initial heads\n", " `Out.Q` : ndarray, shape: (Nt, Nz, Ny, Nx), [L3/T]\n", " net inflow in all cells during time step, inactive cells have 0\n", " Q during time step i, use Out.Q[i]\n", " `Out.Qs` : ndarray, shape: (Nt, Nz, Ny, Nx), [L3/T]\n", " release from storage during time step.\n", " `Out.Qx : ndarray, shape: (Nt, Nz, Ny, Nx-1), [L3/T] \n", " intercell flows in x-direction (parallel to the rows)\n", " `Out.Qy` : ndarray, shape: (Nt, Nz, Ny-1, Nx), [L3/T] \n", " intercell flows in y-direction (parallel to the columns)\n", " `Out.Qz` : ndarray, shape: (Nt, Nz-1, Ny, Nx), [L3/T] \n", " intercell flows in z-direction (vertially upward postitive) \n", " \n", " TO 161024\n", " '''\n", "\n", " import pdb\n", " \n", " # define the named tuple to hold all the output of the model fdm3\n", " Out = namedtuple('Out',['t', 'Phi', 'Q', 'Qs', 'Qx', 'Qy', 'Qz'])\n", " Out.__doc__ = \"\"\"fdm3 output, , containing fields `t`, `Phi`, `Q`, `Qs`, `Qx`, `Qy` and `Qz`\\n \\\n", " Use Out.Phi, Out.Q, Out.Qx, Out.Qy and Out.Qz\n", " or\n", " Out.Phi[i] for the 3D heads of time `i`\n", " Out.Q[i] for the 3D flows of time step `i`\n", " Notice the difference between time and time step\n", " The shape of Phi is (Nt + 1,Nz, Ny, Nx)\n", " The shape of Q, Qs is (Nt, Nz, Ny, Nx)\n", " For the other shapes see docstring of fdm_t \n", " \"\"\" \n", " \n", " if gr.axial:\n", " print('Running in axial mode, y-values are ignored.')\n", "\n", " if kx.shape != gr.shape:\n", " raise AssertionError(\"shape of kx {0} differs from that of model {1}\".format(kx.shape,gr.shape))\n", " if ky.shape != gr.shape:\n", " raise AssertionError(\"shape of ky {0} differs from that of model {1}\".format(ky.shape,gr.shape))\n", " if kz.shape != gr.shape:\n", " raise AssertionError(\"shape of kz {0} differs from that of model {1}\".format(kz.shape,gr.shape))\n", " if Ss.shape != gr.shape:\n", " raise AssertionError(\"shape of Ss {0} differs from that of model {1}\".format(Ss.shape,gr.shape))\n", " \n", " active = (IBOUND>0).reshape(gr.Nod,) # boolean vector denoting the active cells\n", " inact = (IBOUND==0).reshape(gr.Nod,) # boolean vector denoting inacive cells\n", " fxhd = (IBOUND<0).reshape(gr.Nod,) # boolean vector denoting fixed-head cells\n", "\n", " # reshaping shorthands\n", " dx = np.reshape(gr.dx, (1, 1, gr.Nx))\n", " dy = np.reshape(gr.dy, (1, gr.Ny, 1))\n", "\n", " # half cell flow resistances\n", " if not gr.axial:\n", " Rx1 = 0.5 * dx / ( dy * gr.DZ) / kx\n", " Rx2 = Rx1\n", " Ry1 = 0.5 * dy / (gr.DZ * dx) / ky\n", " Rz1 = 0.5 * gr.DZ / ( dx * dy) / kz\n", " else:\n", " # prevent div by zero warning in next line; has not effect because x[0] is not used \n", " x = gr.x.copy(); x[0] = x[0] if x[0]>0 else 0.1* x[1]\n", " \n", " Rx1 = 1 / (2 * np.pi * kx * gr.DZ) * np.log(x[1:] / gr.xm).reshape((1, 1, gr.Nx))\n", " Rx2 = 1 / (2 * np.pi * kx * gr.DZ) * np.log(gr.xm / x[:-1]).reshape((1, 1, gr.Nx))\n", " Ry1 = np.inf * np.ones(gr.shape)\n", " Rz1 = 0.5 * gr.DZ / (np.pi * (gr.x[1:]**2 - gr.x[:-1]**2).reshape((1, 1, gr.Nx)) * kz)\n", " \n", " # set flow resistance in inactive cells to infinite\n", " Rx1[inact.reshape(gr.shape)] = np.inf\n", " Rx2[inact.reshape(gr.shape)] = np.inf\n", " Ry1[inact.reshape(gr.shape)] = np.inf\n", " Ry2 = Ry1\n", " Rz1[inact.reshape(gr.shape)] = np.inf\n", " Rz2 = Rz1\n", " \n", " # conductances between adjacent cells\n", " Cx = 1 / (Rx1[ :, :, :-1] + Rx2[:, : , 1:]) \n", " Cy = 1 / (Ry1[ :, :-1,: ] + Ry2[:, 1:, : ])\n", " Cz = 1 / (Rz1[:-1, :, : ] + Rz2[1:, :, : ])\n", "\n", " # storage term, variable dt not included\n", " Cs = (Ss * gr.Volume / epsilon).ravel()\n", " \n", " # cell number of neighboring cells\n", " IE = gr.NOD[ :, :, 1:] # east neighbor cell numbers\n", " IW = gr.NOD[ :, :, :-1] # west neighbor cell numbers\n", " IN = gr.NOD[ :, :-1, :] # north neighbor cell numbers\n", " IS = gr.NOD[ :, 1:, :] # south neighbor cell numbers\n", " IT = gr.NOD[:-1, :, :] # top neighbor cell numbers\n", " IB = gr.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", " )), (gr.Nod,gr.Nod))\n", "\n", " A = -A + sp.diags(np.array(A.sum(axis=1))[:,0]) # Change sign and add diagonal\n", " \n", " #Initialize output arrays (= memory allocation)\n", " Nt = len(t)-1\n", " Out.Phi = np.zeros((Nt+1, gr.Nod)) # Nt+1 times\n", " Out.Q = np.zeros((Nt , gr.Nod)) # Nt time steps\n", " Out.Qs = np.zeros((Nt , gr.Nod))\n", " Out.Qx = np.zeros((Nt, gr.Nz, gr.Ny, gr.Nx-1))\n", " Out.Qy = np.zeros((Nt, gr.Nz, gr.Ny-1, gr.Nx))\n", " Out.Qz = np.zeros((Nt, gr.Nz-1, gr.Ny, gr.Nx))\n", " \n", " # reshape input arrays to vectors for use in system equation\n", " FQ = R(FQ); HI = R(HI); Cs = R(Cs)\n", " \n", " # initialize heads\n", " Out.Phi[0] = HI\n", " \n", " # solve heads at active locations at t_i+eps*dt_i\n", " \n", " Nt=len(t) # for heads, at all times Phi at t[0] = initial head\n", " Ndt=len(np.diff(t)) # for flows, average within time step\n", " \n", " for idt, dt in enumerate(np.diff(t)):\n", " \n", " it = idt + 1 \n", " \n", " # this A is not complete !!\n", " RHS = FQ - (A + sp.diags(Cs / dt))[:,fxhd].dot(Out.Phi[it-1][fxhd]) # Right-hand side vector\n", "\n", " Out.Phi[it][active] = spsolve( (A + sp.diags(Cs / dt))[active][:,active],\n", " RHS[active] + Cs[active] / dt * Out.Phi[it-1][active])\n", " \n", " # net cell inflow\n", " Out.Q[idt] = A.dot(Out.Phi[it])\n", "\n", " Out.Qs[idt] = -Cs/dt * (Out.Phi[it]-Out.Phi[it-1])\n", "\n", "\n", " #Flows across cell faces\n", " Out.Qx[idt] = -np.diff( Out.Phi[it].reshape(gr.shape), axis=2) * Cx\n", " Out.Qy[idt] = +np.diff( Out.Phi[it].reshape(gr.shape), axis=1) * Cy\n", " Out.Qz[idt] = +np.diff( Out.Phi[it].reshape(gr.shape), axis=0) * Cz\n", "\n", " # update head to end of time step\n", " Out.Phi[it] = Out.Phi[it-1] + (Out.Phi[it] - Out.Phi[it-1]) / epsilon\n", "\n", " # reshape Phi to shape of grid\n", " Out.Phi = Out.Phi.reshape((Nt,) + gr.shape)\n", " Out.Q = Out.Q.reshape( (Ndt,) + gr.shape)\n", " Out.Qs = Out.Qs.reshape((Ndt,) + gr.shape)\n", "\n", " return Out # all outputs in a named tuple for easy access" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Examples\n", "\n", "Here we'll work out a few axially symmetric examples to verify this model using analytical solutions.\n", "\n", "\n", "We will also compute the darwdown in a multi-layer aquifer system.\n", "\n", "We'll keep truly axially symmetric 3D flow for the next chapter, after we introduced the stream function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Preparatory work\n", "\n", "As always we set the path to our own modules and import the required modules and general packages.\n", "\n", "Each time after we edited one or more of our modules, we have to `reload` them. This is why `reload` is imported." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "myModules = './modules/'\n", "\n", "#Adding the path to our modules to the pythonpath\n", "import sys\n", "if not myModules in sys.path:\n", " sys.path.append(myModules)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# import the required general packages\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from importlib import reload\n", "import fdm_t # first time import\n", "reload(fdm_t) # in case we edited fdm_t above we need to reload\n", "\n", "# allow inline plotting\n", "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "def inpoly(...)\n" ] } ], "source": [ "# import our own modules and packages, when they exist\n", "import fdm_t # from current directory\n", "import mfgrid # path has been added to sys.path above\n", "import mfetc\n", "import mfexceptions as err" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "def inpoly(...)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# when we have been editing files, make sure to reload\n", "reload(fdm_t)\n", "reload(mfgrid)\n", "reload(mfetc)\n", "reload(err)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A well in a confined (or unconfined) infinite aquifer (Theis)\n", "\n", "The most famous analytic transient groundwater solution is that of the Theis well, a fully penetrating well in a homogeneous confined aquifer, that starts pumping at a constant rate at t=0. We will now use our model in axially symmetric mode to compute this drawdown numerically and compare it with the analytical solution.\n", "\n", "The analytical solution of the drawdown $s$ Theis well reads\n", "\n", "$$ s = \\frac Q {4 \\pi kD} W(u), \\,\\,\\,\\, u = \\frac {r^2 S} { 4 kD t} $$\n", "\n", "$W(-)$ is called the Theis well function. It is a regular mathematical function known at the `exponential integral`\n", "\n", "$$ W(u) = \\intop_u^\\infty \\frac {e^{-y}} y dy $$\n", "\n", "This exponential integral is available in Python as\n", "\n", "from scipy.special import expi, defined as\n", "\n", "$$ expi(w) = \\intop_{-\\infty}^w \\frac {e^{\\nu}} \\nu d\\nu $$\n", "\n", "To convert the `Well` function `w` to the `expi` function as defined in Python.\n", "\n", "\n", "$$\n", "W(u) =\\intop_{y=u}^{y=\\infty} \\frac {e^{-y}} y dy\n", "= \\intop_{\\nu=-\\infty}^{\\nu=-u} \\frac {e^\\nu} \\nu d\\nu\n", "= \\intop_{y=\\infty}^{y=w} \\frac {e^y} y dy\n", "= expi(w) = expi(-u)\n", "$$\n", "\n", "Which allows us to just use `expi(-u)` for the analytical solution. It may practially be implemented by defining an anonymous function (lambda function or marco) as follows:\n", "\n", "W = lambda u: scipy.linalg.expi(-u)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.0379295765381134" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.special import expi\n", "def W(u): return -expi(-u)\n", "W(0.01) # check if it works" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running in axial mode, y-values are ignored.\n" ] }, { "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" } ], "source": [ "# Extract the heads\n", "plt.figure()\n", "plt.setp(plt.gca(), 'xlabel','x [m]', 'ylabel', 'head [m]' , 'xscale', 'log')\n", "for it in range(len(t)):\n", " plt.plot(gr.xm[gr.xm>0], Out.Phi[it][-1, iy(0), gr.xm>0], 'b.')\n", " plt.plot(gr.xm[gr.xm>0], fi[it][gr.xm>0], 'r')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As can be seen, the numerical and analyticl solutions agree. By changing the variable axial to False or True one can choose between a flat cross sectional model and an axisymmetrical model. The correct axis and analytical solutions will then be shown in the figure. The burden of changing variable for both cases is completely carried by the grid object and the way in which the conductances are computed in fdm3.\n", "\n", "To show that the water balance matches we compute below the total flow into the model." ] }, { "cell_type": "code", "execution_count": 13, "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": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute and show the Hantush type curves together with the Theis type curve\n", "plt.figure()\n", "plt.setp(plt.gca(), 'xscale','log', 'yscale','log',\n", " 'xlabel', '1/u','ylabel','Wh(u)',\n", " 'title', 'Hantush type curves',\n", " 'ylim',(1e-5, 1e2), 'xlim', (0.1, 1e6))\n", "\n", "# Values of r/L\n", "Rho = [0.01, 0.05, 0.1, 0.5, 1., 2., 3.]\n", "\n", "# values for U\n", "U = 1/np.logspace(-1, 6, 71)\n", "\n", "# Theis\n", "plt.plot(1/U, -scipy.special.expi(-U), label='Theis')\n", "# Hantush for each value of r/L\n", "for ir, rho in enumerate(Rho):\n", " plt.plot(1/U, Wh(U,rho), label='r/L={0}'.format(rho)) \n", "plt.legend(fontsize='small')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now try to reproduce these Hantush type curves with our 3D transient finite difference model and compare them with the analytical solution.\n", "\n", "To produce the curves, we make sure that $\\frac Q {4 \\pi kD} = 1$. We also make sure that $1/u$ will vary from 0.1 to $10^6$. And we make sure that we produce curves for the same set of values for $r/\\lambda$." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running in axial mode, y-values are ignored.\n" ] }, { "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 = $('