{ "cells": [ { "cell_type": "markdown", "id": "45d1dfdc", "metadata": {}, "source": [ "# Check Response Functions\n", "This notebooks checks the step response functions by numerically integrating the impulse response functions." ] }, { "cell_type": "code", "execution_count": null, "id": "938b129f-c9ac-418d-879e-f2eba52298d4", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import quad\n", "import pastas as ps" ] }, { "cell_type": "markdown", "id": "9a50691a-85cb-487b-87ad-a238224cdda2", "metadata": {}, "source": [ "### Gamma" ] }, { "cell_type": "code", "execution_count": null, "id": "c5a802df-52d5-49c7-95c0-aa026b94b4f6", "metadata": {}, "outputs": [], "source": [ "ps.Gamma.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "7b8bf54b-44ac-444c-9f30-4fad27b5a546", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "n = 1.5\n", "a = 50\n", "p = [A, n, a]\n", "\n", "gamma = ps.Gamma()\n", "tmax = gamma.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = gamma.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(gamma.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "41d2e835-bf3a-47e7-9f6b-7d0dce2d9c6b", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum, \"--\", label=\"numerical\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] }, { "cell_type": "markdown", "id": "94fb2160-3637-43c7-b8f1-25d508c9b485", "metadata": {}, "source": [ "### Exponential" ] }, { "cell_type": "code", "execution_count": null, "id": "c4508596-cf12-4ed4-b74a-567b25a88ee5", "metadata": {}, "outputs": [], "source": [ "ps.Exponential.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "77ca80bb-c678-46c1-aab1-f0e2ec60ebdf", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "a = 50\n", "p = [A, a]\n", "\n", "exponential = ps.Exponential()\n", "tmax = exponential.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = exponential.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(exponential.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "6172dd00-763c-4f25-b4f0-3c34a15b3c33", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum, \"--\", label=\"numerical\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] }, { "cell_type": "markdown", "id": "5d007aa2-7c99-4c3c-9c2a-f413a8f5a5a3", "metadata": {}, "source": [ "### Hantush" ] }, { "cell_type": "code", "execution_count": null, "id": "750373d0-402a-478c-b1f4-b3457cdb4812", "metadata": {}, "outputs": [], "source": [ "ps.Hantush.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "35790f67-b518-4367-90c3-2d5a425da0c2", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "a = 50\n", "b = 2\n", "p = [A, a, b]\n", "\n", "hantush = ps.Hantush()\n", "tmax = hantush.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = hantush.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(hantush.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "e63a25d6-1b30-41a2-a02d-d4d013925c30", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum, \"--\", label=\"numerical\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] }, { "cell_type": "markdown", "id": "d5d822c8-f954-4e7c-96c6-de24895d4da8", "metadata": {}, "source": [ "### Polder" ] }, { "cell_type": "code", "execution_count": null, "id": "c2dd0ae6-fa59-482b-94b1-74b4b31e0a8a", "metadata": {}, "outputs": [], "source": [ "ps.Polder.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "2eaa17db-8cba-4731-926e-4f8479987063", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "a = 100\n", "b = 0.25\n", "p = [A, a, b]\n", "\n", "polder = ps.Polder()\n", "tmax = polder.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = polder.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(polder.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "419cf91b-9a53-440d-958d-fc085ec0a4b8", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum, \"--\", label=\"numerical\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] }, { "cell_type": "markdown", "id": "cd0facf9-763d-47a6-b722-ebf50b02a2d1", "metadata": {}, "source": [ "### Four-parameter function" ] }, { "cell_type": "code", "execution_count": null, "id": "57d08799-3632-4bb6-84a1-509f72a5861a", "metadata": {}, "outputs": [], "source": [ "ps.FourParam.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "83601f06-0737-4ae9-b7eb-0b113a90bc85", "metadata": {}, "outputs": [], "source": [ "A = 1 # impulse response implemented for A=1 only\n", "n = 1.5\n", "a = 50\n", "b = 10\n", "p = [A, n, a, b]\n", "\n", "fourparam = ps.FourParam(quad=False) # use simple integration\n", "tmax = fourparam.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = fourparam.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(fourparam.impulse, 0, t[i], args=(p))[0]\n", "stepnum = (\n", " stepnum / quad(fourparam.impulse, 0, np.inf, args=p)[0]\n", ") # four param is scaled at the end" ] }, { "cell_type": "code", "execution_count": null, "id": "095f83cb-9a04-4651-8031-7f2f84445c54", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum, \"--\", label=\"numerical\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] }, { "cell_type": "markdown", "id": "d0396351-cd54-4666-b847-862ee442a032", "metadata": {}, "source": [ "### Double exponential function" ] }, { "cell_type": "code", "execution_count": null, "id": "07790a6a-8ec8-4b95-9160-a832b2aceed2", "metadata": {}, "outputs": [], "source": [ "ps.DoubleExponential.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "42e26e95-9afb-4595-b170-0f5ee6e0a186", "metadata": {}, "outputs": [], "source": [ "A = 5 # impulse response implemented for A=1 only\n", "a = 10\n", "b = 50\n", "f = 0.4\n", "p = [A, f, a, b]\n", "\n", "doubexp = ps.DoubleExponential()\n", "tmax = doubexp.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = doubexp.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(doubexp.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "8299b68e-fe97-45cf-8b5b-7f61534d6a60", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum, \"--\", label=\"numerical\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] }, { "cell_type": "markdown", "id": "f5b5fa2c", "metadata": {}, "source": [ "### Kraijenhoff\n", "\n", "#### Kraijenhoff van de Leur\n", "\n", "##### Impulse Response\n", "from [A study of non-steady groundwater flow with special reference to a reservoir-coefficient (1958)](https://edepot.wur.nl/422032) formula 2\n", "\n", "$ \\theta(t) = \\frac{4N}{\\pi S} \\sum_{n=1,3,5...}^\\infty \\left( \\frac{1}{n} \\exp{\\left( {-n^2\\frac{\\pi^2T}{SL^2} t} \\right)} \\sin \\left(\\frac{n\\pi x}{L}\\right) \\right) $\n", "\n", "##### Step Response\n", "\n", "The step response is obtained by taking the integral of the impulse response function\n", "\n", "$ \\Theta(t) = \\frac{4 N}{\\pi S} \\sum_{n=1,3,5...}^\\infty \\frac{1}{n^3} \\left(\\frac{SL^2}{\\pi^2 T} - \\frac{SL^2}{\\pi^2 T} \\exp\\left(-n^2\\frac{\\pi^2T}{SL^2}t\\right)\\right) \\sin \\left(\\frac{n\\pi x}{L}\\right) $\n", "\n", "$ \\Theta(t) = \\frac{4 N L^2}{\\pi^3 T} \\sum_{n=1,3,5...}^\\infty \\frac{1}{n^3} \\left(1 - \\exp\\left(-n^2\\frac{\\pi^2T}{SL^2}t\\right)\\right) \\sin \\left(\\frac{n\\pi x}{L}\\right)$\n", "\n", "And $\\sum_{n=1,3,5...}^\\infty n = \\sum_{n=0}^\\infty (2n+1)$ gives:\n", "\n", "$ \\Theta(t) = \\frac{4 N L^2}{\\pi^3 T} \\sum_{n=0}^\\infty \\frac{1}{(2n+1)^3} \\left(1 - \\exp\\left(-(2n+1)^2\\frac{\\pi^2T}{SL^2}t)\\right)\\right) \\sin \\left(\\frac{(2n+1)\\pi x}{L}\\right)$\n", "\n", "\n", "Kraijenhoff van de Leur takes $\\frac{x}{L}=\\frac{1}{2}$ as the middle of the domain. \n", "\n", "#### Bruggeman \n", "\n", "from Analytical Solutions of Geohydrological Problems (1999) formula 133.15\n", "\n", "##### Step Response\n", "\n", "$ \\Theta(t) = \\frac{-N}{2T}\\left(x^2 - \\frac{1}{4}L^2\\right) - \\frac{4NL^2}{\\pi^3T} \\sum_{n=0}^\\infty\n", "\\frac{(-1)^n}{(2n + 1)^3} \\cos\\left(\\frac{(2n+1)\\pi x}{L}\\right)\n", "\\exp\\left(-\\frac{(2n+1)^2\\pi^2 T}{SL^2}t\\right)\n", "$\n", "\n", "$ \\Theta(t) = \\frac{-NL^2}{2T}\\left(\\left(\\frac{x}{L}\\right)^2 - \\frac{1}{4}\\right) - \\frac{4NL^2}{\\pi^3T} \\sum_{n=0}^\\infty\n", "\\frac{(-1)^n}{(2n + 1)^3} \n", "\\exp\\left(-\\frac{(2n+1)^2\\pi^2 T}{SL^2}t\\right) \\cos\\left(\\frac{(2n+1)\\pi x}{L}\\right)\n", "$\n", "\n", "$ \\Theta(t) = \\frac{-NL^2}{2T}\\left(\\left(\\frac{x}{L}\\right)^2 - \\tfrac{1}{4}\\right) \\left(1 - \\frac{8}{\\pi^3 \\left(\\frac{1}{4} - \\left(\\frac{x}{L}\\right)^2\\right)} \\sum_{n=0}^\\infty\n", "\\frac{(-1)^n}{(2n + 1)^3} \n", "\\exp\\left(-\\frac{(2n+1)^2\\pi^2 T}{SL^2}t\\right) \\cos\\left(\\frac{(2n+1)\\pi x}{L}\\right) \\right)\n", "$\n", "\n", "Note that $x=0$ is the middle of the domain for Bruggeman.\n", "\n", "#### Pastas Implementation\n", "\n", "In Pastas the Bruggeman response function is computed and the parameters are transformed to:\n", "\n", "Scale parameter (such that the gain is always $A$):\n", "\n", "$A = \\frac{-NL^2}{2T}\\left(\\left(\\frac{x}{L}\\right)^2 - \\tfrac{1}{4}\\right)$\n", "\n", "Reservoir coefficient (also known as $j$ in Kraijenhoff):\n", "\n", "$a = \\frac{SL^2}{\\pi^2 T}$\n", "\n", "Location in the domain:\n", "\n", "$b = \\frac{x}{L}$\n", "\n", "\n", "Such that the step response becomes:\n", "\n", "$ \\Theta(t) = A\\left(1 - \\frac{8}{\\pi^3(\\frac{1}{4} - b^2)} \\sum_{n=0}^\\infty \\frac{(-1)^n}{(2n+1)^3} \\cos\\left((2n+1)\\pi b\\right)\\exp\\left(-\\frac{(2n+1)^2t}{a}\\right) \\right)$\n", "\n", "Taking the derivative gives the impulse response:" ] }, { "cell_type": "code", "execution_count": null, "id": "798ab1b1", "metadata": {}, "outputs": [], "source": [ "ps.Kraijenhoff.impulse" ] }, { "cell_type": "code", "execution_count": null, "id": "c3835b9e", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "a = 10\n", "b = 0.25\n", "p = [A, a, b]\n", "\n", "khoff = ps.Kraijenhoff()\n", "tmax = khoff.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = khoff.step(p)\n", "stepnum_brug = np.zeros(len(t))\n", "\n", "for i in range(1, len(t)):\n", " stepnum_brug[i] = quad(khoff.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "e2c870cd", "metadata": {}, "outputs": [], "source": [ "plt.plot(t[1:], step, label=\"analytic\")\n", "plt.plot(t, stepnum_brug, \"--\", label=\"numerical bruggeman\")\n", "plt.xlabel(\"time (d)\")\n", "plt.ylabel(\"step (m)\")\n", "plt.grid()\n", "_ = plt.legend() # try to show figure in readthedocs" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "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.9.12" }, "vscode": { "interpreter": { "hash": "29475f8be425919747d373d827cb41e481e140756dd3c75aa328bf3399a0138e" } } }, "nbformat": 4, "nbformat_minor": 5 }