{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test polygon area sink"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[1, 2], z=[10, 5, 4, 0], c=2)\n",
    "xy = [(-50, 0), (50, 0), (50, 80), (-50, 80)]\n",
    "p1 = tfs.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy,\n",
    "    kaq=[1, 2],\n",
    "    z=[10, 5, 4, 0],\n",
    "    c=[2],\n",
    "    topboundary=\"conf\",\n",
    "    N=0.01,\n",
    "    order=5,\n",
    "    ndeg=3,\n",
    ")\n",
    "rf = tfs.Constant(ml, xr=0, yr=-1000, hr=2)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(win=[-100, 100, -100, 100], ngr=100, layers=[0], levels=20, labels=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-200, 200, 100)\n",
    "h = ml.headalongline(x, 3)\n",
    "plt.plot(x, h[0])\n",
    "plt.plot(x, h[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checks for numerical derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# recharge inside polygon (should be 0.01)\n",
    "x = 20\n",
    "y = 60\n",
    "d = 0.01\n",
    "d2hdx2 = (ml.head(x + d, y) - 2 * ml.head(x, y) + ml.head(x - d, y)) / (d**2)\n",
    "d2hdy2 = (ml.head(x, y + d) - 2 * ml.head(x, y) + ml.head(x, y - d)) / (d**2)\n",
    "d2hdx2 + d2hdy2\n",
    "aqin = ml.aq.inhomlist[0]\n",
    "print(\"recharge from numerical derivative: \", np.sum(aqin.T * (d2hdx2 + d2hdy2)))\n",
    "h = ml.head(x, y)\n",
    "print(\"leakage from aq0 to aq1 from head difference: \", (h[1] - h[0]) / aqin.c[1])\n",
    "print(\n",
    "    \"leakage from aq0 to aq1 from num. derivative: \",\n",
    "    aqin.T[1] * (d2hdx2[1] + d2hdy2[1]),\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}