{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](\"AW&H2015.tiff\")
\n",
"
\n",
"
\n",
"
\n",
"# Problem P3.4 Theis solution with two pumping\n",
"\n",
"In this notebook, we will work through one tutorial based the Theis solution for transient pumping, and investigate the superposition of drawdown from two interfering pumping wells. Two wells fully penetrate a 20-m-thick confined aquifer that is isotropic and homogeneous (Fig. P3.1). Storativity is estimated to be 2 x 10-5. The hydraulic conductivity is 100 m/d. The confining unit is composed of very low permeability material and is approximated as impermeable. Both wells have a radius of 0.5 m and are pumped continuously at a constant rate for 30 days; well A is pumped at 4000 L/min and well B is pumped at 12,000 L/min. Before pumping, the head is 100 m everywhere in the problem domain. The 800 m by 500 m problem domain in Fig. P3.1 is the near-field region of a problem domain that extends over many tens of square kilometers so that the aquifer effectively is of infinite extent and the composite cone of depression does not reach the boundaries after 30 days of pumping.\n",
"\n",
"We simpflied it to look like this:\n",
"\n",
"![](\"P3.4_figure.tiff\")
\n",
"\n",
"Below is an iPython Notebook that builds a Theis function and plots results. \n",
"\n",
"[Acknowledgements: This tutorial was created by Randy Hunt and all failings are mine. The exercise here is modeled after example iPython Notebooks developed by Chris Langevin and Joe Hughes for the USGS Spring 2015 Python Training course GW1774]
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"'/Users/rjhunt1/GitHub/Chapter_3_problems-1'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Problem 3.4, page 107 Anderson, Woessner and Hunt (2015)\n",
"\n",
"# import Python libraries/functionality for use in this notebook\n",
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.special\n",
"import sys, os\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"\n",
"# return current working directory\n",
"os.getcwd()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Name of model path: /Users/rjhunt1/GitHub/Chapter_3_problems-1/P3-4_Theis\n",
"Creating model working directory.\n"
]
}
],
"source": [
"# Set the name of the path to the model working directory\n",
"dirname = \"P3-4_Theis\"\n",
"datapath = os.getcwd()\n",
"modelpath = os.path.join(datapath, dirname)\n",
"print 'Name of model path: ', modelpath\n",
"\n",
"# Now let's check if this directory exists. If not, then we will create it.\n",
"if os.path.exists(modelpath):\n",
" print 'Model working directory already exists.'\n",
"else:\n",
" print 'Creating model working directory.'\n",
" os.mkdir(modelpath)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"'/Users/rjhunt1/GitHub/Chapter_3_problems-1/P3-4_Theis'"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"os.chdir(modelpath)\n",
"os.getcwd()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Define an function, class, and object for Theis Well analysis\n",
"\n",
"def well_function(u):\n",
" return scipy.special.exp1(u)\n",
"\n",
"def theis(Q, T, S, r, t):\n",
" u = r ** 2 * S / 4. / T / t\n",
" s = Q / 4. / np.pi / T * well_function(u)\n",
" return s\n",
"\n",
"class Well(object):\n",
" def __init__(self, x, y, rate, name):\n",
" self.x = float(x)\n",
" self.y = float(y)\n",
" self.rate = rate\n",
" self.name = name\n",
" self.swell = None\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Parameters needed to solve Theis\n",
"r = 500 # m\n",
"T = 2000 # m^2/d (100 m/d Kh x 20 m thick)\n",
"S = 0.00002 # unitless\n",
"t = 30. # days\n",
"\n",
"#Q = pumping rate # m^3/d - but we'll enter it below in the well info"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Well information\n",
"well_list =[]\n",
"well_obj = Well(250, 250, 5760, \"Well A\") # 4000 L/min = 5760 m^3/d\n",
"well_list.append(well_obj)\n",
"well_list.append(Well(550, 250, 17280, \"Well B\")) # 12000 L/min = 17280 m^3/d\n",
"\n",
"# Grid information as requested in problem\n",
"x = np.linspace(0, 800., 50) # x-direction 0 to 800 m, 50 m increments\n",
"y = np.linspace(0, 500., 50) # y-direction 0 to 500 m, 50 m increments\n",
"xgrid, ygrid = np.meshgrid(x, y) # make a grid with these coordinates"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### We want to explore drawdown as a function of time\n",
"So, set up an array of times to evaluate, and loop over them. Also, we can specify a distance from each well at which to calculate the curve of drawdown over time."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.\n",
" 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29.\n",
" 30.]\n"
]
}
],
"source": [
"times = np.linspace(0.,30.,31) # linear interpolation of time from 0 to 30 days, make 30 increments days at 0.5\n",
"rdist = 25 # this sets the distance to plot drawdown over time\n",
"print times"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### We will want to normalize our plots\n",
"Let's figure out the maximum drawdown to use for setting our colorbar on the plots."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"15.4068299862\n"
]
}
],
"source": [
"#let's find the maximum drawdown\n",
"drawdown_grid_max = np.zeros(xgrid.shape, dtype=float)\n",
"for well_obj in well_list:\n",
" r = ((well_obj.x - xgrid)**2 + (well_obj.y - ygrid) ** 2) ** 0.5\n",
" s_max = theis(well_obj.rate, T, S, r, times[-1])\n",
" drawdown_grid_max += s_max\n",
"max_drawdown = np.max(drawdown_grid_max)\n",
"print max_drawdown"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Clobber the PNG output files\n",
"We will make a sequence of PNG files from which to render the movie. It's a good idea to delete the existing ones from our folder. \n",
"\n",
" ** Do this kind of thing with caution!** \n",
"\n",
" ** You are deleting files and will _not_ be asked to confirm anything!**"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for cf in os.listdir(os.getcwd()):\n",
" if cf.endswith('.png'):\n",
" os.remove(cf)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Loop over time and make figures\n",
"We will make a figure with the drawdown contours over the whole grid. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"working on time 30.0: 100.00% complete"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD3CAYAAAA0Vx7KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXFWd9/HPlwCGkJAIYYBAIJFhSUAhjERkGRsHIS6A\nOiLBQREchhnWB32UxDU8jriMjuiADqPARGQVhUlE2ZRAECSgCYYkDSIECEtYEzbjBPJ7/rinkkql\nqrqqu7rqVtX3/XrVq+9ybt1fd93761Pn3nuOIgIzM8uvjVodgJmZVedEbWaWc07UZmY550RtZpZz\nTtRmZjnnRG1mlnNO1GZmfZA0RVKvpD9KOqtCme+m9fdKmlTPtn1xojYzq0LSEOA8YAowEThG0oSS\nMu8B/joidgH+Cfh+rdvWwonazKy6ycCDEbE0IlYDVwBHlpQ5ApgJEBF3AaMkbVvjtn1yojYzq257\n4LGi+WVpWS1lxtSwbZ+cqM3Mqqu1nw0NVgAbD9Ybm5nlhaS6OjWKiOKk+zgwtmh+LFnNmCpldkhl\nNqlh2z45UZtZl7ilxnIHly64B9hF0jjgCeBo4JiSMrOAU4ErJO0HrIiI5ZKeq2HbPjlRm5lVERGv\nSToVuAEYAlwYEUsknZTWXxARv5D0HkkPAq8Ax1fbtt4Y5G5OzazTZU0ftdeoS5o+Ws4XE83Mcs6J\n2sws55yozcxyzonazCznnKjNzHLOidrMLOecqM3Mcs6J2sws55yozcxyzom6y0laI+lNrY6jUSTt\nKOklSbl6ssxsIJyom0zSUkmvSnpR0guSfiPpJCeWDUmaIemSPsoslfTOwnxEPBoRIyLHfSNI+rGk\nJ9Mx8JCkz5Ws/7s0dNMrkn4tacdWxWr54ETdfAG8LyK2AHYEvgacBVxYaQNJ/pwqCwaxH+BB8lVg\nfDoG3g2cJmkKgKTRwE+BzwFvJOu57cpWBWr54ATQQhHxUkTMJuv68DhJEwEk/bek70v6haSXgR5J\n75U0X9JKSY9K+lLhfSTNlPTJNL19as44Oc3vnLpaLJT9tKQnJC2TdEJxPJJGSvqRpKdTTfVzhZq+\npEck7ZOm/yHtY0Ka/4Ska9L0DElXpZhelHSfpL+p9DeQ9J30+6yUdI+kA9PyKcB04OjUlDG/zLaX\nkP2zm53K/F9J41JsG6UycyR9OX1zeUnSLEmjJV2a9jlP0k5F77m7pJskPZdqtUfV/IHWKCIWRcSq\nokWvAU+n6Q8C90XETyPif4EZwF6Sdm10HNY+nKhzICLuJutM/KCixccAX46I4cBvgJeBYyNiJPBe\n4F8kFcZemwP0pOl3AA8Bf1s0fxusTX6fAg4Bdk0/i/0HMAIYn7b7GKm7xjL7+FP6WZifU/Q+hwOX\nAyPJ+uk9r8qvPw/Yi6z2eBnwE0mbRsT1wDnAFakpY1LphhHxUeBRsm8oIyLimxX2cTRwLNkQSDsD\nd5J9g9kSWAJ8CUDS5sBNwI+BrYGpwPcqDUYq6Xup+arca0GV37mw7SvAIuBfI+L3adUewL1Fv+Or\nwIPAntXezzqbE3V+PEGWOAqujYg7ASLiLxFxa0QsSvMLyQbJLCTK24ADU+33IOAbwAFp3TuAW9P0\nh4GLImJxSgDFtfIhZAltekS8EhGPAN8CPpqK3Fq0vwPJvr4X5v+2aB8AcyPi+tRO/GOyRFxWRFwa\nES9ExJqI+HfgDcBuhbAYeLNGABdHxMMR8SLwS+CBiPh1RLwO/AQo/BN4H/BwRMxM8SwAfgaUrVVH\nxMkR8cYKr72rBhVxMjCc7J/lv0qanFZtDrxYUvzFVNa6lBN1fuwAPJ+mg/UHxETS2yTdkpolVgAn\nAVsBRMSfyDor35ssUf8ceCJ9XS5OotuVvO+jRdOjyYYNeqRkfWEgztuAg5SNrDyELMEdkJoNRqak\nVrC8aPpVYGildvbUXLFY0gpJL5DVwkeXKzsAxfGsYl0zQ2G+kAR3At5WXDMGPgJs0+B4AIjMHLK/\nZWHUj5eBLUqKjgReGowYrD04UeeApH3JRiu+vUqxy4BrgR0iYhTwn6z/+d1KVvPbJCKeSPMfJ2tS\nKCTRJ8nadAuKp58FVgPjStYvA4iIB8mS7mnArRHxEvAU8E/A3KJtar7bQtJBwKeBoyJiVES8EVjJ\nulp0Le9V790d1co/Sva7FdeMR0TEKeUKS/rP1O5d7rWwjpg2IftHC1lTyNpvIKk5Zue03FpA0kWS\nllf6TFNlY356LZT0mqRRtWxbKyfq1ihcoNtC0vvI2nMvKTRtUP7r/nDghYj43/Q1+SOsn3RuJRuz\n7bY0PyfNzy26Ve0q4OOSJkgaRlHTR2oGuAr4iqThqaZ8JlnTRek+CjX0OSXzlWKvZATZhbRnJW0q\n6YusX5t8ChhXuKBZwXKyRFaNKkyXug7YVdKxkjZJr30l7V6ucET8c0rk5V5vLhuItLWkqZI2lzRE\n0mFk/2D/JxW5BthT0gclDSX7jBZExAN9/I42eC4GplRaGRHfjIhJ6TrKdGBORKyoZdtaOVG3xmxJ\nL5LV4KaTtQUfX7Q+2LDmdzLw/9J2X2DDW7ZuI0vmhUT9G2CzonnSBbpzgV8DDwC/KtnPaWQ1u4fI\nasmXkh1oBbeW7KN0vlLslWqx16fXA8BS4M+s3xzzk/TzOUn3VHiPrwKfT00Vn6ywvyiZLrs+fUs4\nlOwi4uNk30C+CmxaYd/9EcA/k31TeQ74MvDRdEGZiHgW+HvgK2RNYW9N8ViLRMRc4IUai3+ErOLV\nn20rqmnMRElLyS5ovA6sjojJkrYkSxY7kZ1kHy78F5E0HTghlT89Im4caKBmZv2lAY6ZqGwU8dmV\nvimlMsPIrgHtXFSjrmnbvtQ6CnkAPRHxfNGyacBNEfENSWel+WnK7gU+GphIdiHqZkm7RsSa/gZp\nZjZgE3rKL39lDrw6Z938s/3ew+HA7cVJulFqTdSwYdveEay7PWsmWXvlNOBI4PKIWA0sVTZ8+mTg\ntwML1cxsEGzek70Knj27v+80laJmj0aqtY06yGrG90g6MS3bJiIKtz0tZ90tTGNIdwoky1h3i5eZ\nWceRNJLsVtj/6atsf9Raoz4gIp6UtDVwk6Te4pUREVkbUEXrreujrJnZekrbjJtJ0uVkrQejJT1G\ndifOJimuC1Kx9wM3RMSfK2y7Vdr2ixFRfIG+JjUl6oh4Mv18RlmfDpOB5ZK2jYinJG3HuocIHgfG\nFm2+Q1q2niWxU+mipjtvxgpOnTGq62NoRBy7z36k70J9mHEZzPjIgN+mI+IYSAy9hzfm3MrLsTlB\nAz+2BiIijqmhzEyyJuC6t61Fn00fkoZJGpGmNye7fWkhWR8Ox6Vix5E9jEFaPjXdFzse2IWsPwfr\nQLvPfqQhSdoax59J56mlRr0NcE165mBj4NKIuDHd13qVpE+Qbs8DiIjFkq4CFpM9zHBynvsGtvo5\nCbSHwufUqBq2tU6fiToiHibrQ6J0+fNs2PtaYd05ZD2f5drknqGtDiEXMUBtcQx2gu7p912mjZWH\nOBoZw+6zH+lXss7LsWk1PvDS8J1KkYc2aquNa9Cdox1r1xP0yIAvJkoKJtSY65aopRcvy6nnPmrr\nMk7QncfNIe3JfX3YBnwxqvP5820vrlHbWj55u4tr1+3DNWoDnKS7mT/7/HONusv5JDVw7TrvXKPu\nUm6HtnJ8TOSTE3WXcYK2vvj4yB8n6i7iE9Bq5WMlX9xG3QV80ll/uN06P1yj7nBO0jZQPoZaz4m6\nQ7kt2hrJx1JrOVF3IJ9UNhh8XLWOE3UHcS3aBls3Hl+Shkq6S9ICSYslfbVK2X0lvSbpg2l+rKRb\nJC2SdJ+k0/sTgxN1h+jGE8hao9uOtYhYBRwcEXsDbwEOlnRgaTlJQ4CvA9ezbjDw1cCZEbEHsB9w\niqQJ9cbgRN3mXIu2Vui2Yy4iXk2TmwJDgOfLFDsNuBp4pmi7pyJiQZp+GVhCNgB4XZyo21i3nSyW\nL910/EnaSNICYDlwS0QsLlm/PXAk8P20aIPOryWNAyYBd9W7f99H3aa66SSx/Orv6DEt8X8rLO+d\nA/fPWTe/ZMMiEbEG2FvSSOAGST0RUbQR5wLTIiKUjVu43sADkoaT1bbPSDXruniElzbjBG15NJjJ\numEjvFxYY677RPURXiR9AfhzRHyzaNlDrEvOo4FXgRMjYpakTYCfA7+MiHP7E7+bPtqIk7TlVScf\nm5JGSxqVpjcD3gXMLy4TEW+KiPERMZ6s5vwvKUkLuBBY3N8kDU7UbaOTTwTrDB18jG4H/Dq1Ud8F\nzI6IX0k6SdJJfWx7AHAs2Z0i89NrSr0BuI26DXTwCWCWexGxENinzPILKpQ/vmj6dhpQIXaNOuec\npK2d+HgdHE7UOeX7o61d+bhtPCfqHPKBbu3Ox3Bj1ZSoJQ1JjeCz0/yWkm6S9ICkGwtXRNO66ZL+\nKKlX0qGDFXin8gFuncLHcuPUWqM+A1jMuqdtpgE3RcSuwK/SPJImAkcDE4EpwPckudZeIx/YZlZO\nn0lU0g7Ae4Afsu6G7iOAmWl6JvD+NH0kcHlErI6IpcCDwORGBtypnKStE/m4boxaarvfBj4NrCla\ntk1ELE/Ty4Ft0vQYYFlRuWXA9gMNstP5YLZO5uN74KreRy3pfcDTETFfUk+5MunZ9mrPZpZdd96M\nFWunJ/cMZXLP0L6j7UA+iM3WN2/OKubNWdXqMHKlrwde9geOkPQeYCiwhaRLgOWSto2IpyRtBzyd\nyj8OjC3afoe0bAOnzhhVbnFXcZK2blFP502lFbfzz145WGG1japNHxHx2YgYm55fnwr8OiI+CswC\njkvFjgOuTdOzgKmSNpU0HtgFmDc4obc3J2nrNj7m+6/eR8gLzRhfA66S9AlgKfBhgIhYLOkqsjtE\nXgNOjlZ0z5dzPmDNrB7u5rTJnKSt29XbJWreujltBd/j3ERO0mbWH07UZtZUrrDUz4m6SXxwmll/\nOVE3gZO02fp8TtTHiXqQ+YA0a3+Slkr6Q+qcboNbjiXtLulOSaskfarM+vU6tquXE/UgcpI2q6zN\nzo8AeiJiUkSU67/oOeA04Jtl1sGGHdvVxYnazKw2FW/Zi4hnIuIeYPUGG5Xv2K4uTtSDpM1qC2Yt\n0UbnSQA3S7pH0ol1bluuY7u6eHDbQdBGB59Z19jphN6yy1fNmceqOeuanSv0LHJARDwpaWvgJkm9\nETG3r33W0rFdLZyozayl6umwaTAM7ZnM0J51zc4rzz5/gzIR8WT6+Yyka8j62e8zUVO+Y7sfRcTH\n6onRTR8N5tq0WWeRNEzSiDS9OXAosLBS8eKZCh3b1ZWkwYm6oZykzfon5+fONsBcSQuAu4CfR8SN\nkk6SdBKApG0lPQacCXxe0qOShpd5r37d9eGmDzOzKiLiYWDvMssvKJp+ivX74i/3PrcCt/YnBteo\nGyTnNQIza2NO1GaWC67sVOZE3QA+wMxsMLmNeoCcpLvIhndt1e6UhkVhXciJ2qyagSTnSu/jpF1R\nq++pzisn6gFwbbpDNSo51/L+TtpWAydqs4LBTtDV9umEbVX4YqIZtCZJ52n/OeJvqhtyjbqffDB1\niDwlSNeurQLXqK175SlJF8trXNYyTtT94Np0B8h7Msx7fNZUVRO1pKGS7pK0QNJiSV9Ny7eUdJOk\nByTdKGlU0TbTJf1RUq+kQwf7FzCry/m0TxJslzgHgStD66uaqCNiFXBwROwNvAU4WNKBwDTgpojY\nFfhVmkfSROBoYCIwBfieJNfaLR/aMfG1Y8zWcH0m0Yh4NU1uCgwBXgCOAGam5TOB96fpI4HLI2J1\nRCwFHiTrYLtj+D99m2rnhNfOsVtD9JmoJW2U+mFdDtwSEYuAbSJieSqynKy/VoAxwLKizZcB2zcw\nXrPu5GTd1fq8PS8i1gB7SxoJ3CDp4JL1IalaZ9hl1503Y8Xa6ck9Q5ncM7S2iM3q1SlJ7ny64ta9\nOQuz17O/W9F34SZJ1+F+COxBltNOiIjfFq0fDfwY2JYsr34zIv47rVsKvAi8DqyOiLpbGWq+jzoi\nVkq6DvgbYLmkbSPiKUnbAU+nYo+zfufZO6RlGzh1xqhyi3PNzR5tqFOSdBfpeXP2gpX0Hr4T559d\nYbjZ5voO8IuI+JCkjYHNS9afCsyPiOkpad8v6ccR8RpZYu+JiOf7u/O+7voYXbijQ9JmwLuA+cAs\n4LhU7Djg2jQ9C5gqaVNJ44FdgHmYtUInJulO/J1yLrUmHBQRFwFExGsRUfrf40lgizS9BfBcStJr\n32YgMfRVo94OmJnu3NgIuCQifiVpPnCVpE8AS4EPp19gsaSrgMXAa8DJEdGvMcLMrIIuaQLJkfHA\nM5IuBvYCfgecUXSjBcAPgF9LegIYQcqJSQA3S3oduCAiflBvAFUTdUQsBPYps/x54JAK25wDnFNv\nIGYN5ZqnlfgYl9RU7ssbLtqYLA+eGhF3SzqX7JbkLxaV+SywICJ6JO0M3CRpr4h4CTggIp6UtHVa\n3hsRc+uJ3X191Mjt05YrrlU3zNI5j/DInKrn9zJgWUTcneavJj07UmR/4CsAEfEnSQ8DuwH3RMST\nafkzkq4hu2XZidq6nGvTVodxPTsxrmfdYAW3nX37euvTTROPSdo1Ih4ga01YVPI2vWn5byRtQ5ak\nH5I0DBgSES9J2hw4FDi73hidqM3M+nYacKmkTYE/ASdIOgkgIi4ga+69WNK9ZNfzPhMRz0t6E/Az\nSZDl20sj4sZ6d+5EbTYI5twAPYcN8k66oPkjL02OEXEvsG/J4guK1j8LHF5mu4eAvQe6fydq6yxN\nbvaYc0Pf6wY9YVvHc4dJHeDMH8J3Zq2bP+xLcOJ56+Y/dSF8+38qb//xc+Gnd2TTPZ+F3z1Yvtyz\nL8ImH4ALrh94zJ2gUpI+n+xqU8G+N8CJ962b/1QvfHtp5ff9+EL46VPZdM88+F2Z5z165sHut8Ok\nM2DiKfCDKv8wrP05UXeAAyfCHb3Z9Jo18NxLsPjRdevvvB8OmFB5e2nd3fhS9irnJ7fDlH3g8tsa\nEnZbq1aT3pN1V5rWkD07fOfj67a5cyUcUOXBXLHuMyieLi1z2Vtg/nfgN1+Hs2bCa6/X+UtY23Ci\nrkFe2skqeftucGdK1IsehT13hBGbwYqX4S+rYcljsM/OWU2557Pw1k/ClC/BUy/Ut58r5sK/HgtP\nr4THn2v879EuqiVpyDqDWJyml5I9LTEMeBm48QZY8jLss0VWU+6ZB2+9E6bcA0/9pb44IoDz4cVX\nYfhQGOKzuWO5jboDjNkKNh4Cjz2TJey3754l0jvvhy02g7eMy8qd9l8w+/Ow1RZw5Vz43CVw4em1\n7eOxZ7IEvdd4+ND+2faffH/f2zVVTm7LG03WH/DTZDXricCzaXoY8JYRWbnTemH2JNhqU7jySfjc\nH+HCPWvbRwD/sBDeIPjjLfCdEyt/E7L250TdIfbfPWv+uKMXPnkkPP483LEERm6eNXvc/3hW2z7k\nC1n519fAmC1rf/8rb88SNMBRB8AJ381hom6CvmrTBXsA95El56NYl6g3B3Z4Ae5/FRa9DIfck5V/\nPWDMG2qPo9D0sc8W8OxHYf/PwGH7wI5b1/4e1j6cqDvEARPgN0tg4SPw5nEwdmv45jUwchic8K7s\na/IeO8Id3+jf+19+GyxfAT+ek80/+QI8+AT89ZhG/QadZU+yRP0Q8Cbgr4ArgeHAu0mfx3C4420D\n39foLbKmrbvud6LuVG7V6hD7T4Cf3w1bjci+Ar9xOKx4JWv+2H932HUMPLMSfpvasle/tv4Fx2oe\neBxeWQXLLoaHf5i9pv29LypWswfwW2AkWe13BFkb9aK0btfN4Zn/hd+mLpdXr4HFL9e3j0J3Z6/+\nBeY/BH+9XWNit/xxjbpD7LljdrfHsT3rlr1lXHYSb5naRK+eBqf/F6x8NbtD4MwjYOKOfb/3FXPh\ng29ff9nf7w9T/w2+MLVRv0FnGU92t0dxz2U7A38h6wNz043g6r3g9F5Y+Rq8FnDmTjBxeO37+IeF\nsNlG8JfFcPzfwaSdG/kbWJ6oFb2QSoolsVPfBXMi73d9WNKEi4m1tlH3paEPwXT404k6AiJiQJdK\nJcUX4rM1lf2yzhnw/hrNTR9mZjnnRG1Wh0bUhP1IudXLidrMLOecqGvQe3j7tKd3tSa11Q6kRuza\ntPWHE7VZP/Qn4Q5Kku7wC4mWcaI266d6Eq9r0jYQTtRmA9BzWPUk3Nd6aw+SpktaJGmhpMskvaFk\nfY+klZLmp9fni9aNknS1pCWSFkvar979+4EX6yyn0JLOmVqSjN3s0RSSxgEnAhMi4i+SrgSmAjNL\nit4aEUeUeYvvAL+IiA9J2pisy5e6uEZtZrmVkwv5LwKrgWEp0Q4DHi9TboOHZCSNBA6KiIsAIuK1\niCgzFER1TtTWeVzTtAaKiOeBbwGPAk8AKyLi5tJiwP6S7pX0C0kT0/LxwDOSLpb0e0k/SCOT18WP\nkNfBj5K3kZz0TT1ouuSfUe/hOzFBjzTkEfJac07p/iTtDMwGDgJWAj8Bro6IS4vKjABej4hXJb0b\n+E5E7CrprcCdwP4Rcbekc4EXI+KL9cTfZxu1pLHAj8h6agzgvyLiu5K2JOu5cSeygSw+HBEr0jbT\ngROA14HT+zM8utmAtKit2trPvDmrmDdnVbUibwXuiIjnACT9DNgfWJuoI+KloulfSvpeypHLgGUR\ncXdafTUwrd4Ya7mYuBo4MyIWSBoO/E7STcDxwE0R8Q1JZ6WdT0tV/qPJBrbYHrhZ0q4Rsabe4Mys\njC6pTTfL5J6hTO4Zunb+/LM3aELuBb4gaTNgFVmniPOKC0jaBng6IkLSZLLWiufTusdSDnwgbbuI\nOvXZRh0RT0XEgjT9MrCELAEfwbqrnjOBwngfRwKXR8TqiFgKPAhMrjcwswFzQrMGiIh7yVoV7gH+\nkBb/QNJJkk5K8x8CFkpaAJxLdldIwWnApZLuBd4CnFNvDHXdnpduU5kE3AVsExHL06rlwDZpegxZ\nn+kFy8gSe9vrPXwnt1Nba/mfT0tExDeA0vGRLihafz4VGttSot93IPuvOVGnZo+fAmdExEsqGkkz\nVferXZXcYN15M1asnS796mHWMJ3UVt0lSXrOwuz17G4j4Xcr+t6gC9SUqCVtQpakL4mIa9Pi5ZK2\njYinJG1HNugyZPcXji3afAfK3HN46oxR/Y/arB6dkKy7JEkD9Lw5e/UenuWIMm3GXafPNmplVecL\ngcURcW7RqlnAcWn6OODaouVTJW0qaTywCyUN72ZN10WJzjpPLQ+8HAAcCxxc9Bz7FOBrwLskPQC8\nM80TEYuBq4DFwC+Bk6MVN2sPkpw8KWX90a7Jul3jtobps+kjIm6nckI/pNzCiDiHflzZNBt07dYM\n4iRt+BHyfnGtus21S/JrlzgHgc+x9TlRW3fKcxI8hXzHZ03nRG3dK48JMW/xWC44UfeTv5p1kLwk\n7DzEkAM+tzbkgQPMClpxodHJ2WrgRD0AfqS8AxUnzsFM2k7QVgcnarNKGp20nZytn5yoB8i16i5R\nKclWSuBOyv3i9unynKjNBsIJ2ZrAd300gGsBZgPn86gyJ2ozs5xzom4Q1wbMOpOkiyQtl7SwaNmW\nkm6S9ICkGyVt0G+zpLGSbpG0SNJ9kk4vWndFUSd3D0uaXy0GJ+oGcrI265+cnzsXA1NKlk0jGzN2\nV+BXlB+wtjDe7B7AfsApkiYARMTUiJgUEZPI+vr/abUAnKjNzKqIiLnACyWLK40ZW7xdufFmxxSX\nSf39fxi4vFoMTtQNlvOagZk1RqUxY8sqGW+22EHA8oj4U7XtfXveIPC91Wa1a1blptI5WRijsb/6\nGjM2jTd7Ndl4sy+XrD4GuKyvfThRm1lXK4zRWHD2FTVtVmnM2PUUjTf746LxZgvrNgY+AOzT187c\n9DFI3ARi1rc2Pk8qjRm7VpXxZgsOAZZExBN97cyJehC18UFoZomky4E7gN0kPSbpeCqMGStpjKTr\n0qblxpt9d9FbH00fFxHXxtCKcWclxZLoniTm9mqzDdVakZmgR4gIDWRfkiJm1Vj2CAa8v0ZzjdrM\nLOecqJvATSBm6/M5UR8n6ibxgWlm/eVE3URO1mY+D/qjz0Rdb4ckkqZL+qOkXkmHDlbg7coHqXUz\nH//9U0uNuuYOSSRNJLvlZGLa5nuSXGsv4YPVzOrRZxKts0OSI4HLI2J1RCwFHgQmNybUzuJkbd3G\nx3z/9be2W6lDkjHAsqJyy4Dt+7mPjucD18xqMeC+PvrqkAQou+68GSvWTk/uGcrknqEDDaUtuQMn\n6wb1VErmzVnFvDmrBjGa9tPfRF2pQ5LHgbFF5XZIyzZw6owNBkToWk7W1snq/eZYWnE7/+yVjQ6p\n7fS36aNShySzgKmSNpU0HtgFmDewELuDm0GsE/m4boxabs+ruUOSiFgMXAUsBn4JnByt6EykTfmg\nNrNy3ClTDrkZxDpBoyoe7pTJTybmkmvW1u58DDeWE3VO+UA3swIn6hxzsrZ25OO28Zyoc84HvbWT\nTjxeK/R3dJSkRZJel1RxzENJZ0haKOk+SWcULZ8saV4a9eVuSftWi8GJug30Hr5TR54A1lk6+Bgt\n19/RQrKBaW+rtJGkPYF/BPYF9gLeJ2nntPobwBciYhLwxTRfkRN1G+ngE8HaXCcfm+X6O4qI3oh4\noI9NdwfuiohVEfE6cCvwwbTuSWBkmh5FhQcDC5yo20wnnxDWnnxMVnQfcFDqFnoY8F6yp7Uh63H0\nW5IeBf4NmF7tjQbc14c1nx85t7xoqyR9fvnFc57PXo0WEb2Svg7cCLwCzAdeT6svBE6PiGskHQVc\nBLyr0nv5gZc254RtrdKsJN2wB14Oq7HsDRs+8CJpHDA7It5csvwW4FMR8fsaYjgHeDQi/lPSixGx\nRVouYEVEjKy0rZs+2lxb1WisY/i4W0/FfyKS/ir93JHs4uNladWDkt6Rpt8JVG3vdtNHByicNK5d\nWzN0W5JO/R29Axgt6THgS8DzwH8Ao4HrJM2PiHdLGgP8ICLemza/WtJWwGqyvo9eTMv/CThf0huA\nP6f5yjFiXAg8AAAIU0lEQVS46aOzOFnbYGpFks5D00eruemjw/ieaxssPq5ax4m6Q/mkskby8dRa\nTtQdzLVrawQfQ63nRN0FfKJZf/nYyQcn6i7h2rXVy8dLfjhRdxmffNYX/1PPHyfqLuQT0SrxcZFP\nTtRdzAnbivlYyC8/mWh+srHLOUHnn2vUtpZP2O7jz7w9uEZt63Htujs4QbcX16itLLdfdy5/ru1n\nUGrUkqYA5wJDgB9GxNcHYz82+FzD7hxO0O2r4TVqSUOA88gGg5wIHCNpQqP30wjz5qxqdQi5iAH6\njqNQwx7Mk33Owr7LNEMe4mhkDP393PJybNrgNH1MBh6MiKURsRq4AjhyEPYzYHk4EPMQA9QXx2Al\n7DwkSMhHHI2KYSCfU16OTRucpo/tgceK5pcBbxuE/ViLuVkkv9zM0VkGo0bd/JEIrKV84TE//Fk0\nnqSLJC2XtLBo2Zcl3StpgaRfSRpbYdtRkq6WtETSYkn7peV7SbpT0h8kzZI0omoMjR7hJQUyIyKm\npPnpwJriC4qSnMzNrGatHOFF0kHAy8CPCoPbShoRES+l6dOAvSLiH8vsdyZwa0RcJGljYPOIWCnp\nbuCTETFX0vHA+Ij4YqWYBqPp4x5glzRq7xPA0cAxxQXyNsyNmVklKZmOK1n2UtHscODZ0u0kjQQO\niojj0javASvT6l0iYm6avhm4HqiYqBve9JGCORW4AVgMXBkRSxq9HzOzVpL0FUmPAscBXytTZDzw\njKSLJf1e0g8kDUvrFkkq3GRxFFC26WTtvloxuK2ZWTNJilsqrFuQXgUz2fBbf6pRzy40fZSsmwbs\nFhHHlyx/K3AnsH9E3C3pXODFiPiipN2A7wJbAbOA0yNidKX4/Qi5mXW1vdOrYGb9b3EZ8Isyy5cB\nyyLi7jR/NTANICLuBw4DkLQr8N5qO2j6I+SSpkjqlfRHSWcN4n7KXandUtJNkh6QdKOkUUXrpqeY\neiUd2qAYxkq6RdIiSfdJOr1FcQyVdFe6Qr1Y0ldbEUd63yGS5kua3cIYlqar7fMlzWtFHGXuBnhb\nC2LYLf0NCq+Vkk5vQRzT0zmyUNJlkt7QiuOizph3KZo9EphfWiYingIeS4kY4BBgUdp+6/RzI+Dz\nwPer7jAimvYie6T8QWAcsAnZN44Jg7Svg4BJwMKiZd8APpOmzwK+lqYnplg2SbE9CGzUgBi2BfZO\n08OB+4EJzY4jvfew9HNj4LfAgS2K45PApcCsVnwm6b0fBrYsWdbsY2MmcELRZzKyFX+Long2Ap4k\nayttWhzpfR4C3pDmryRr821oDEDcUuMLiJJtLye7MeJ/yZ4ROYGsdrwwxfJT4K9S2THAdUXb7gXc\nDdwL/AwYmZafTpYP7gfO6TP+Rn7YNfyx3g5cXzQ/DZg2iPsbx/qJuhfYJk1vC/Sm6enAWUXlrgf2\nG4R4riX7r9qyOIBh6cDZo9lxADuQXeE+mKy9ryWfCVmi3qpkWdPiIEvKD5VZ3srj4lBgbgv+Flum\nZPVGsn9Ys4F3NToGBpCo8/BqdtNHuacWt2/i/reJiOVpejmwTZoek2IZtLjSxYhJwF2tiEPSRpIW\npP3dEhGLWhDHt4FPA2uKlrXiMwngZkn3SDqxBXGUuxtg8ybHUGoqWc2RZsYREc8D3wIeJau1roiI\nm5oZQztodqLOzS0mkf2brRZPw2KVNJzs69EZsf79l02LIyLWRMTeZLXav5V0cDPjkPQ+4OmImA+U\nvY++iZ/JARExCXg3cIqyBxqaGcfGwD7A9yJiH+AV0kWmJsawlqRNgcOBn2ywk8E/LnYG/g/Zt98x\nwHBJxzYzhnbQ7ET9OOvfLziW9f87DrblkrYFkLQd8HSFuHZIywZM0iZkSfqSiLi2VXEURMRK4Drg\nb5ocx/7AEZIeJqu5vVPSJU2OAYCIeDL9fAa4hqwjsWbGUe5ugH2Ap1p0XLwb+F36e0Bz/xZvBe6I\niOciewbjZ2RNpK36W+RSsxP12qcW03/xo8nuIWyWWWQXKkg/ry1aPlXSppLGA7sA8wa6M0kCLgQW\nR8S5LYxjdOGquaTNyNoA5zczjoj4bESMjYjxZF+zfx0RH21mDACShin1q5CaGw4luyjUzL9FpbsB\nZjcrhhLHsK7Zo7C/ZsXRC+wnabN0vhxC9qBcq/4W+dTsRnGy/973k12tnT6I+ym9Uns82YWLm4EH\ngBuBUUXlP5ti6gUOa1AMB5K1xy4gS4zzyfrpbnYcbwZ+n+L4A/DptLypcRS99ztYd9dHs/8W41n3\njMN9hWOwBXFscDdAKz4PYHOyx59HFC1r9t/iM2T/qBaS3Q2zSaNjoM0vJvrJRDPreKryZGKpg8lf\nf0QeM9HMLOecqM3Mcs6J2sws55yozcxyzonazCznnKjNzHLOidrMLOecqM3Mcs6J2sws55yozcyq\nUPnRomZIWlY0Os6UKtuvN6pRWnZUGtXmdUn79BWDE7WZWXUXk/XRUyyAf4+ISel1fZXtzyDraKq4\nv46FwAeA22oJwInazKyKiJgLvFBmVZ/9gUjaAXgP8MPi8hHRGxEP1BqDE7WZWf+cJuleSRcWD75b\notyoRnVzojYzq9/3ybrM3ZtsUOBvlRaoZVSjWm08kI3NzNrFwX0XqVlEFEacQdIPyQY6KFUY1eg9\nwFBgC0k/ioiP1bs/16jNrONFhOp59fV+aXiwgg+QXRws3We5UY3KJek+9+dEbWZWhaTLgTuA3SQ9\nJukE4OuS/iDpXrIRi85MZcdIuq7CW62960PSByQ9BuwHXCfpl1Vj8AgvZmb55hq1mVnOOVGbmeWc\nE7WZWc45UZuZ5ZwTtZlZzjlRm5nlnBO1mVnOOVGbmeXc/wdkBPcMaO7joQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Note that this section of code is saving figures for animation - not plotting them!\n",
"\n",
"from IPython.display import clear_output\n",
"# to make our plots of drawdown over time a one point, we can\n",
"# predefine the response as np.nan. That way, when we plot incrementally\n",
"# as we calculate through time, only the times for which calculations\n",
"# have been made will appear using plt.plot()\n",
"for well_obj in well_list:\n",
" well_obj.swell = np.ones_like(times)*np.nan\n",
" \n",
"\n",
"# using \"enumerate\" we get both the iterant (t) and a counter (i) \n",
"for i,t in enumerate(times):\n",
" # the following stuff just writes out a status message to the screen\n",
" clear_output()\n",
" perc_done = (i/float(len(times)-1)) * 100\n",
" sys.stdout.write('working on time {0}: {1:2.2f}% complete'.format(t,\n",
" perc_done))\n",
" if i < len(times):\n",
" sys.stdout.flush()\n",
" # here's the end of the silly shenanigans of plotting out status to the screen \n",
"\n",
" # now we calculate the drawdown for each time.\n",
" drawdown_grid = np.zeros(xgrid.shape, dtype=float)\n",
" for well_obj in well_list:\n",
" r = ((well_obj.x - xgrid)**2 + (well_obj.y - ygrid) ** 2) ** 0.5\n",
" s = theis(well_obj.rate, T, S, r, t)\n",
" well_obj.swell[i] = (theis(well_obj.rate, T, S, rdist, t))\n",
" drawdown_grid += s\n",
" \n",
" \n",
" # drawdown contour map (map view)\n",
" plt.subplot(1, 3, 1, aspect='equal')\n",
" im = plt.contourf(xgrid, \n",
" ygrid, \n",
" drawdown_grid, \n",
" np.linspace(0,max_drawdown,10))\n",
" # optional color bar configuration\n",
" divider = make_axes_locatable(plt.gca())\n",
" cax = divider.append_axes(\"right\", \"5%\", pad=\"3%\")\n",
" plt.colorbar(im, cax=cax).ax.invert_yaxis()\n",
" for well_obj in well_list:\n",
" plt.text(well_obj.x, well_obj.y, well_obj.name)\n",
" plt.title('Drawdown at time = {0:.0f}'.format(t))\n",
" \n",
"\n",
" \n",
"# Let's finish with a drawdown only plot --> make a second set of figures with only the \n",
" # make a plot \n",
" plt.subplot(1, 1, 1, aspect='equal')\n",
" im = plt.contourf(xgrid, \n",
" ygrid, \n",
" drawdown_grid, \n",
" np.linspace(0,max_drawdown,10))\n",
" plt.colorbar().ax.invert_yaxis()\n",
" for well_obj in well_list:\n",
" plt.text(well_obj.x, well_obj.y, well_obj.name)\n",
" plt.title('Drawdown at time = {0:.0f}'.format(t))\n",
" plt.savefig('s_only{0}.png'.format(i))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Let's make an animation!"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ffmpeg_path is: ffmpeg\n"
]
}
],
"source": [
"# for execution robustness, we need to determine where ffmpeg lives\n",
"# in general, you probably won't need to bother\n",
"import platform\n",
"from subprocess import check_output\n",
"if 'Windows' in platform.platform():\n",
" if '64bit' in platform.architecture()[0]:\n",
" ffmpeg_path = os.path.join(binpath, 'ffmpeg.exe')\n",
" else:\n",
" ffmpeg_path = os.path.join(binpath, 'win32', 'ffmpeg.exe')\n",
"else:\n",
" #Assume it is in path on macos\n",
" ffmpeg_path = 'ffmpeg'\n",
"print 'ffmpeg_path is: ', ffmpeg_path\n",
"\n",
"figfiles = ['s_only%d.png']\n",
"anmfiles = ['Theis_movie1.mp4']\n",
"\n",
"# note the tricky way we can iterate over the elements of \n",
"# two lists in pairs using zip (if you wanted to add more plots)\n",
"for figfile,anmfile in zip(figfiles,anmfiles):\n",
" try:\n",
" os.remove(anmfile)\n",
" print 'Deleted the existing animation: ', anmfile\n",
" except:\n",
" pass\n",
" # now we do a system call, making the movie using command line arguments\n",
" # for ffmpeg\n",
" output = check_output([ffmpeg_path, \n",
" '-f', 'image2', \n",
" '-i', figfile, \n",
" '-vcodec', 'libx264',\n",
" '-pix_fmt', 'yuv420p',\n",
" anmfile])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Finally, we can embed the drawdown movie into the notebook\n",
"We could also just look at it outside of the notebook by looking in the working subdirectory P3-4_Theis. Note too that you can manually move through time by grabbing the time slider in the playbar. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"