{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b2c5b32e",
   "metadata": {},
   "source": [
    "# Test wells with resistance\n",
    "\n",
    "This notebook checks the implementation for wells with specified screen resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8816e402-3e11-4d2b-8aee-0dd66770caf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31dbe5e3-0ff2-4566-9c3e-51ce6ef4df84",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters single layer\n",
    "k = 10  # m/day\n",
    "c = 1000.0  # resistance top leaky layer in days\n",
    "z = [1.0, 0, -20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef7efed3-ad39-4acf-8444-36b6f02285e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res)\n",
    "ml.solve()\n",
    "x = 100\n",
    "y = 0\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is {ml.head(x, y).item():.5f}\")\n",
    "print(f\"head inside well: {w.headinside().item():.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acb74ed9-f2a8-4b9b-a488-07852237a134",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0.01\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res)\n",
    "ml.solve()\n",
    "hin = w.headinside()\n",
    "hout = ml.head(w.xc, w.yc, w.layers)\n",
    "Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res\n",
    "x = 100\n",
    "y = 0\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is {ml.head(x, y).item():.5f}\")\n",
    "print(f\"head inside well: {w.headinside().item():.5f}\")\n",
    "print(\"discharge: \", w.discharge())\n",
    "print(\"discharge from equation\", Qcheck)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "468a46e7-7649-47dc-96f2-a1c9145bcf34",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters multiple layers\n",
    "k = [10, 20, 40]  # m/day\n",
    "c = [1000, 500, 800]  # resistance top leaky layer in days\n",
    "z = [1.0, 0, -10, -12, -18, -20, -32]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4032603-f371-44de-b1d1-e64d0e5dd251",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res, layers=[1, 2])\n",
    "ml.solve()\n",
    "x = 100\n",
    "y = 0\n",
    "layer = 2\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is\", ml.head(x, y))\n",
    "print(\"head inside well is\", w.headinside())\n",
    "print(\"discharge:\", w.discharge())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd4d2039-cc87-45dc-bc99-cf4fbc5e1e7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0.01\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res, layers=[1, 2])\n",
    "ml.solve()\n",
    "hin = w.headinside()\n",
    "hout = ml.head(w.xc, w.yc, w.layers)\n",
    "Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res\n",
    "x = 100\n",
    "y = 0\n",
    "layer = 2\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is\", ml.head(x, y))\n",
    "print(\"head inside well is\", w.headinside())\n",
    "print(\"discharge:\", w.discharge())\n",
    "print(\"discharge from equation\", Qcheck)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3ab83f1-1460-40b1-98b4-4f40f8f99d0f",
   "metadata": {},
   "source": [
    "### HeadWell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8847f1ea-385b-437c-91af-9609e9d934a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters\n",
    "k = 10  # m/day\n",
    "c = 1000.0  # resistance top leaky layer in days\n",
    "z = [1.0, 0, -20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b915f1c-284d-4036-9764-7ae4c2105d9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res)\n",
    "ml.solve()\n",
    "x = 100\n",
    "y = 0\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is {ml.head(x, y).item():.5f}\")\n",
    "print(f\"head inside well: {w.headinside().item():.5f}\")\n",
    "print(\"discharge:\", w.discharge())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1135ac1b-b227-43d2-b739-6c23d807574f",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0.01\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res)\n",
    "ml.solve()\n",
    "hin = w.headinside()\n",
    "hout = ml.head(w.xc, w.yc, w.layers)\n",
    "Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res\n",
    "x = 100\n",
    "y = 0\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is {ml.head(x, y).item():.5f}\")\n",
    "print(f\"head inside well: {w.headinside().item():.5f}\")\n",
    "print(\"discharge:\", w.discharge())\n",
    "print(\"discharge from equation\", Qcheck)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f38bf622-01cd-4c80-8b73-9523c728e7e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters multiple layers\n",
    "k = [10, 20, 40]  # m/day\n",
    "c = [1000, 500, 800]  # resistance top leaky layer in days\n",
    "z = [1.0, 0, -10, -12, -18, -20, -32]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f3fe856-8562-46d4-b9e8-6e2cf036317a",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res, layers=[1, 2])\n",
    "ml.solve()\n",
    "x = 100\n",
    "y = 0\n",
    "layer = 2\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is\", ml.head(x, y))\n",
    "print(\"head inside well is\", w.headinside())\n",
    "print(\"discharge:\", w.discharge())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20a96b3d-e7a4-42b6-bd4d-cef40cba65ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0.01\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res, layers=[1, 2])\n",
    "ml.solve()\n",
    "hin = w.headinside()\n",
    "hout = ml.head(w.xc, w.yc, w.layers)\n",
    "Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res\n",
    "x = 100\n",
    "y = 0\n",
    "layer = 2\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is\", ml.head(x, y))\n",
    "print(\"head inside well is\", w.headinside())\n",
    "print(\"discharge:\", w.discharge())\n",
    "print(\"discharge from equation\", Qcheck)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1883e95-af5c-44e8-8f53-e70b3b2fad02",
   "metadata": {},
   "source": [
    "### TargetHeadWell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0a00053-0cfa-4546-b609-bc8e4fdae6e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters single layer\n",
    "k = 10  # m/day\n",
    "c = 1000.0  # resistance top leaky layer in days\n",
    "z = [1.0, 0, -20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6f2430-b1b5-41c9-82d5-836cf86b74a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.TargetHeadWell(ml, 0, 0, rw=0.1, res=res, hcp=-0.4, xcp=100, ycp=200)\n",
    "ml.solve()\n",
    "x = 100\n",
    "y = 200\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is {ml.head(x, y).item():.5f}\")\n",
    "print(f\"head inside well: {w.headinside().item():.5f}\")\n",
    "print(\"discharge:\", w.discharge())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60170120-abbd-43a4-8380-d2b884986c45",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0.01\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.TargetHeadWell(ml, 0, 0, rw=0.1, res=res, hcp=-0.4, xcp=100, ycp=200)\n",
    "ml.solve()\n",
    "hin = w.headinside()\n",
    "hout = ml.head(w.xc, w.yc, w.layers)\n",
    "Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res\n",
    "x = 100\n",
    "y = 200\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is {ml.head(x, y).item():.5f}\")\n",
    "print(f\"head inside well: {w.headinside().item():.5f}\")\n",
    "print(\"discharge:\", w.discharge())\n",
    "print(\"discharge from equation\", Qcheck)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e75c9188-e1de-4dc3-b183-4c797bcf55ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters multiple layers\n",
    "k = [10, 20, 40]  # m/day\n",
    "c = [1000, 500, 800]  # resistance top leaky layer in days\n",
    "z = [1.0, 0, -10, -12, -18, -20, -32]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce1bff47-0ca2-459c-a1a7-d7b156db7c33",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.TargetHeadWell(\n",
    "    ml, 0, 0, rw=0.1, res=res, layers=[1, 2], hcp=-0.4, xcp=100, ycp=200, lcp=1\n",
    ")\n",
    "ml.solve()\n",
    "x = 100\n",
    "y = 200\n",
    "layer = 2\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is\", ml.head(x, y))\n",
    "print(\"head inside well is\", w.headinside())\n",
    "print(\"discharge:\", w.discharge())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7528d7d-88ea-498b-b318-9d598ab08917",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0.01\n",
    "ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary=\"semi\", hstar=0)\n",
    "w = tfs.TargetHeadWell(\n",
    "    ml, 0, 0, rw=0.1, res=res, layers=[1, 2], hcp=-0.4, xcp=100, ycp=200, lcp=1\n",
    ")\n",
    "ml.solve()\n",
    "hin = w.headinside()\n",
    "hout = ml.head(w.xc, w.yc, w.layers)\n",
    "Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res\n",
    "x = 100\n",
    "y = 200\n",
    "layer = 2\n",
    "print(f\"resistance: res={res:.5f}\")\n",
    "print(f\"head at x={x}, y={y} is\", ml.head(x, y))\n",
    "print(\"head inside well is\", w.headinside())\n",
    "print(\"discharge:\", w.discharge())\n",
    "print(\"discharge from equation\", Qcheck)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}