{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synthetic Pumping Test - Calibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from timflow import transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (6, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use observation times from Oude Korendijk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drawdown = np.loadtxt(\"../04pumpingtests/data/oudekorendijk_h30.dat\")\n",
    "\n",
    "tobs = drawdown[:, 0] / 60 / 24\n",
    "robs = 30\n",
    "Q = 788"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=60, z=(-18, -25), Saq=1e-4, tmin=1e-5, tmax=1)\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.1, tsandQ=[(0, 788)], layers=0)\n",
    "ml.solve()\n",
    "rnd = np.random.default_rng(2)\n",
    "hobs = ml.head(robs, 0, tobs)[0] + 0.05 * rnd.random(len(tobs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### See if `timflow.transient` can find aquifer parameters back"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq\", layers=0, initial=100)\n",
    "cal.set_parameter(name=\"Saq\", layers=0, initial=1e-3)\n",
    "cal.series(name=\"obs1\", x=robs, y=0, layer=0, t=tobs, h=hobs)\n",
    "cal.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"rmse:\", cal.rmse())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm = ml.head(robs, 0, tobs, 0)\n",
    "plt.semilogx(tobs, hobs, \".k\")\n",
    "plt.semilogx(tobs, hm[0], \"r\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"correlation matrix\")\n",
    "print(cal.fitresult.covar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit with `scipy.least_squares` (not recommended)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq\", layers=0, initial=100)\n",
    "cal.set_parameter(name=\"Saq\", layers=0, initial=1e-3)\n",
    "cal.series(name=\"obs1\", x=robs, y=0, layer=0, t=tobs, h=hobs)\n",
    "cal.fit_least_squares(report=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibrate parameters in multiple layers\n",
    "Example showing how parameters can be optimized when multiple layers share the same parameter value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=[10.0, 10.0],\n",
    "    z=(-10, -16, -18, -25),\n",
    "    c=[10.0],\n",
    "    Saq=[0.1, 1e-4],\n",
    "    tmin=1e-5,\n",
    "    tmax=1,\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.1, tsandQ=[(0, 788)], layers=1)\n",
    "ml.solve()\n",
    "hobs0 = ml.head(robs, 0, tobs, layers=[0])[0]\n",
    "hobs1 = ml.head(robs, 0, tobs, layers=[1])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.set_parameter(\n",
    "    name=\"kaq\", layers=[0, 1], initial=20.0, pmin=0.0, pmax=30.0\n",
    ")  # layers 0 and 1 have the same k-value\n",
    "cal.set_parameter(name=\"Saq\", layers=0, initial=1e-3, pmin=1e-5, pmax=0.2)\n",
    "cal.set_parameter(name=\"Saq\", layers=1, initial=1e-3, pmin=1e-5, pmax=0.2)\n",
    "cal.set_parameter(name=\"c\", layers=1, initial=1.0, pmin=0.1, pmax=200.0)\n",
    "cal.series(name=\"obs0\", x=robs, y=0, layer=0, t=tobs, h=hobs0)\n",
    "cal.series(name=\"obs1\", x=robs, y=0, layer=1, t=tobs, h=hobs1)\n",
    "cal.fit(report=False)\n",
    "display(cal.parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogx(tobs, hobs0, \".C0\", label=\"obs layer 0\")\n",
    "plt.semilogx(tobs, hobs1, \".C1\", label=\"obs layer 1\")\n",
    "\n",
    "hm = ml.head(robs, 0, tobs)\n",
    "plt.semilogx(tobs, hm[0], \"C0\", label=\"modelled head layer 0\")\n",
    "plt.semilogx(tobs, hm[1], \"C1\", label=\"modelled head layer 1\")\n",
    "\n",
    "plt.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate data for head measured in well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tobs2 = np.hstack((tobs, np.arange(0.61, 1, 0.01)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=60, z=(-18, -25), Saq=1e-4, tmin=1e-5, tmax=1)\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.3, res=0.02, tsandQ=[(0, 788), (0.6, 0)], layers=0)\n",
    "ml.solve()\n",
    "rnd = np.random.default_rng(2)\n",
    "hobs2 = w.headinside(tobs2)[0] + 0.05 * rnd.random(len(tobs2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq\", layers=0, initial=100)\n",
    "cal.set_parameter(name=\"Saq\", layers=0, initial=1e-3)\n",
    "cal.set_parameter_by_reference(name=\"res\", parameter=w.res[:], initial=0.05)\n",
    "cal.seriesinwell(name=\"obs1\", element=w, t=tobs2, h=hobs2)\n",
    "cal.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm = w.headinside(tobs2)\n",
    "plt.semilogx(tobs2, hobs2, \".k\")\n",
    "plt.semilogx(tobs2, hm[0], \"r\");"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}