{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Theis Storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import quad, solve_ivp\n",
    "from scipy.special import exp1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = -100\n",
    "npor = 0.3\n",
    "k = 10\n",
    "H = 10\n",
    "T = k * H\n",
    "Ss = 1e-4\n",
    "S = Ss * H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def headtheis(r, t, T=T, S=S, Q=Q):\n",
    "    u = r**2 * S / (4 * T * t)\n",
    "    h = -Q / (4 * np.pi * T) * exp1(u)\n",
    "    return h\n",
    "\n",
    "\n",
    "def voltheis(r, t, T, S, Q):\n",
    "    u = r**2 * S / (4 * T * t)\n",
    "    h = -Q / (4 * np.pi * T) * exp1(u)\n",
    "    vol = h * 2 * np.pi * r\n",
    "    return vol * S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demonstrate headtheis works correctly\n",
    "quad(voltheis, a=1e-3, b=10000, args=(10, T, S, Q))  # gives 1000 m^3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def volume(r, t, T=T, S=S, Q=Q):\n",
    "    return quad(voltheis, a=1e-3, b=r, args=(10, T, S, Q))[0] + npor * np.pi * (r**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vxytheis(t, xy):\n",
    "    x, y = xy\n",
    "    r = np.sqrt(x**2 + y**2)\n",
    "    u = S * r**2 / (4 * T * t)\n",
    "    Qr = -Q / (2 * np.pi) / r * np.exp(-u)\n",
    "    vr = Qr / (H * npor)\n",
    "    vx = vr * x / r\n",
    "    vy = vr * y / r\n",
    "    return np.array([vx, vy])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vxytheisnew(t, xy):\n",
    "    x, y = xy\n",
    "    r = np.sqrt(x**2 + y**2)\n",
    "    u = S * r**2 / (4 * T * t)\n",
    "    Qr = -Q / (2 * np.pi) / r * np.exp(-u)\n",
    "    vr = Qr / (H * (npor + S * headtheis(x, t)))\n",
    "    vx = vr * x / r\n",
    "    vy = vr * y / r\n",
    "    return np.array([vx, vy])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trace pathline for 10 days\n",
    "path0 = solve_ivp(vxytheis, (1e-5, 10), y0=[1e-5, 0], method=\"DOP853\")\n",
    "R0 = path0.y[0, -1]\n",
    "R0, np.pi * R0**2 * npor * H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path0 = solve_ivp(vxytheisnew, (1e-5, 10), y0=[1e-5, 0], method=\"DOP853\")\n",
    "R1 = path0.y[0, -1]\n",
    "R1, np.pi * R1**2 * npor * H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "headtheis(5, 10) * S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npor * H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = path0.y[0, -1]\n",
    "print(\"R, volume\", R, np.pi * R**2 * npor * H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "quad(headtheis, 1e-5, R, args=(100, T, S, Q))[0] + np.pi * R**2 * npor * H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = 0\n",
    "t1 = 10\n",
    "\n",
    "\n",
    "def vxytheis2(t, xy):\n",
    "    x, y = xy\n",
    "    r = np.sqrt(x**2 + y**2)\n",
    "    u1 = S * r**2 / (4 * T * (t - t0))\n",
    "    u2 = S * r**2 / (4 * T * (t - t1))\n",
    "    Qr = -Q / (2 * np.pi) / r * np.exp(-u1)\n",
    "    if t >= t1:\n",
    "        Qr += 2 * Q / (2 * np.pi) / r * np.exp(-u2)\n",
    "    vr = Qr / (H * npor)\n",
    "    vx = vr * x / r\n",
    "    vy = vr * y / r\n",
    "    return np.array([vx, vy])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path1 = solve_ivp(vxytheis2, (10 + 1e-6, 20), y0=[R, 0], method=\"DOP853\")\n",
    "path1.y[0, -1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vxytheis2(12 + 1e-6, (10, 0))"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}