freyberg_mf6_api_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4115a470",
   "metadata": {},
   "source": [
    "In this notebook, we will see a very simple demonstration of using the MODFLOW-6 API capability for, you guessed, the Freyberg model...\n",
    "\n",
    "First import as usual..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aadb009",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import platform\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import flopy\n",
    "import modflowapi\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4ef5e6e",
   "metadata": {},
   "source": [
    "We can use pre-compiled mf6 binaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b65da487",
   "metadata": {},
   "outputs": [],
   "source": [
    "if \"linux\" in platform.platform().lower():\n",
    "    lib_path = os.path.join(\"bin\",\"linux\",\"libmf6.so\")\n",
    "    exe_path = os.path.join(\"bin\",\"linux\",\"mf6\")\n",
    "    lib_name = os.path.split(lib_path)[-1]\n",
    "    exe_name = \"./\"+os.path.split(exe_path)[-1]\n",
    "elif \"darwin\" in platform.platform().lower() or \"macos\" in platform.platform().lower() :\n",
    "    lib_path = os.path.join(\"bin\",\"mac\",\"libmf6.so\")\n",
    "    exe_path = os.path.join(\"bin\",\"mac\",\"mf6\")\n",
    "    lib_name = os.path.split(lib_path)[-1]\n",
    "    exe_name = \"./\"+os.path.split(exe_path)[-1]\n",
    "else:\n",
    "    lib_path = os.path.join(\"bin\",\"win\",\"libmf6.dll\")\n",
    "    exe_path = os.path.join(\"bin\",\"win\",\"mf6.exe\")  \n",
    "    lib_name = os.path.split(lib_path)[-1]\n",
    "    exe_name = os.path.split(exe_path)[-1]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60c57428",
   "metadata": {},
   "source": [
    "Copy the original model files to a working dir and get the binaries into that dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de5d5e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "org_dir = \"monthly_model_files_org\"\n",
    "work_dir = \"temp\"\n",
    "if os.path.exists(work_dir):\n",
    "    shutil.rmtree(work_dir)\n",
    "shutil.copytree(org_dir,work_dir)\n",
    "shutil.copy2(lib_path,os.path.join(work_dir,os.path.split(lib_path)[-1]))\n",
    "shutil.copy2(exe_path,os.path.join(work_dir,os.path.split(exe_path)[-1]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "900442a9",
   "metadata": {},
   "source": [
    "Now just run standard mf6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9472ac9",
   "metadata": {},
   "outputs": [],
   "source": [
    "c_d = os.getcwd()\n",
    "os.chdir(work_dir)\n",
    "os.system(exe_name)\n",
    "os.chdir(c_d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95a13fa5",
   "metadata": {},
   "source": [
    "Now lets see how the API is used - first just replicate the standard MF6 solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a883ee1",
   "metadata": {},
   "outputs": [],
   "source": [
    "api_dir1 = \"api1\"\n",
    "if os.path.exists(api_dir1):\n",
    "    shutil.rmtree(api_dir1)\n",
    "shutil.copytree(org_dir,api_dir1)\n",
    "shutil.copy2(lib_path,os.path.join(api_dir1,os.path.split(lib_path)[-1]))\n",
    "shutil.copy2(exe_path,os.path.join(api_dir1,os.path.split(exe_path)[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0840b3d9",
   "metadata": {},
   "source": [
    "Create a `ModflowApi` instance for our model and initialize it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "988fb97a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gwf = modflowapi.ModflowApi(os.path.join(api_dir1, lib_name), working_directory=api_dir1)\n",
    "gwf.initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "576b5356",
   "metadata": {},
   "source": [
    "Now march thru the stress periods and solve as usual..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e24980fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get current sim time\n",
    "ctime = gwf.get_current_time()\n",
    "# get ending sim time\n",
    "etime = gwf.get_end_time()\n",
    "# max number of iterations\n",
    "max_iter = gwf.get_value(gwf.get_var_address(\"MXITER\", \"SLN_1\"))\n",
    "# let's do it!\n",
    "while ctime < etime:\n",
    "    # the length of this sim time\n",
    "    dt = gwf.get_time_step()\n",
    "    # prep the current time step\n",
    "    gwf.prepare_time_step(dt)\n",
    "    kiter = 0\n",
    "    # prep to solve\n",
    "    gwf.prepare_solve(1)\n",
    "    # the current one-based stress period number\n",
    "    stress_period = gwf.get_value(gwf.get_var_address(\"KPER\", \"TDIS\"))[0]\n",
    "    time_step = gwf.get_value(gwf.get_var_address(\"KSTP\", \"TDIS\"))[0]\n",
    "    # solve until converged\n",
    "    while kiter < max_iter:\n",
    "        if gwf.solve(1):\n",
    "            print(\"flow stress period,time step {0},{1} converged with {2} iters\".\\\n",
    "                  format(stress_period, time_step, kiter))\n",
    "            break\n",
    "\n",
    "        kiter += 1\n",
    "    try:\n",
    "        gwf.finalize_solve(1)\n",
    "    except:\n",
    "        pass\n",
    "\n",
    "    gwf.finalize_time_step()\n",
    "    # update current sim time\n",
    "    ctime = gwf.get_current_time()\n",
    "gwf.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f2d7e3",
   "metadata": {},
   "source": [
    "Now let's do something more exciting.  We can use the API interface to build in operational rules for the pumping wells.  If the sfr downstream flow for the terminal reach is too low, we need to pump less, if the flow is high, we can pump more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7335f972",
   "metadata": {},
   "outputs": [],
   "source": [
    "api_dir2 = \"api2\"\n",
    "if os.path.exists(api_dir2):\n",
    "    shutil.rmtree(api_dir2)\n",
    "shutil.copytree(org_dir,api_dir2)\n",
    "shutil.copy2(lib_path,os.path.join(api_dir2,os.path.split(lib_path)[-1]))\n",
    "shutil.copy2(exe_path,os.path.join(api_dir2,os.path.split(exe_path)[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a24fbea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "gwf = modflowapi.ModflowApi(os.path.join(api_dir2, lib_name), working_directory=api_dir2)\n",
    "gwf.initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd2c3392",
   "metadata": {},
   "source": [
    "Now we use march thru stess periods as before, but within the outer iteration loop, we will attempt some \"management\" to link the groundwater extraction rates to surface-water flow.  If the surface-water flow at the terminal downstream reach is less than 2,500, reduce the pumping rates.  Otherwise, if the surface-water flow is greater than 3,000, lets pump more!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd42d51d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# get current sim time\n",
    "ctime = gwf.get_current_time()\n",
    "# get ending sim time\n",
    "etime = gwf.get_end_time()\n",
    "# max number of iterations\n",
    "max_iter = gwf.get_value(gwf.get_var_address(\"MXITER\", \"SLN_1\"))\n",
    "# let's do it!\n",
    "while ctime < etime:\n",
    "    # the length of this sim time\n",
    "    dt = gwf.get_time_step()\n",
    "    # prep the current time step\n",
    "    gwf.prepare_time_step(dt)\n",
    "    kiter = 0\n",
    "    # prep to solve\n",
    "    gwf.prepare_solve(1)\n",
    "    # the current one-based stress period number\n",
    "    stress_period = gwf.get_value(gwf.get_var_address(\"KPER\", \"TDIS\"))[0]\n",
    "    time_step = gwf.get_value(gwf.get_var_address(\"KSTP\", \"TDIS\"))[0]\n",
    "    # solve until converged\n",
    "    while kiter < max_iter:\n",
    "        \n",
    "        # get the sfr dsflow rate for the last reach\n",
    "        addr = [\"DSFLOW\", \"FREYBERG6\", \"SFR_1\"]\n",
    "        wbaddr = gwf.get_var_address(*addr)\n",
    "        sfr_dsflow = gwf.get_value(wbaddr)\n",
    "        \n",
    "        #get the wel extraction rates\n",
    "        \"FREYBERG6/WEL-1                   BOUND\"\n",
    "        addr = [\"BOUND\", \"FREYBERG6\", \"WEl-1\"]\n",
    "        wbaddr = gwf.get_var_address(*addr)\n",
    "        wel_bound = gwf.get_value_ptr(wbaddr)\n",
    "        org = np.sum(wel_bound)\n",
    "        \n",
    "        # if the flow is too low, turn down the wells\n",
    "        fac = 1\n",
    "        if sfr_dsflow[-1] < 2500:\n",
    "            fac = sfr_dsflow[-1]/2500\n",
    "            wel_bound[:] *= fac\n",
    "        \n",
    "        #otherewise if there is extra flow, turn up the wells\n",
    "        elif sfr_dsflow[-1] > 3000:\n",
    "            fac = sfr_dsflow[-1]/3000\n",
    "            wel_bound[:] *= fac\n",
    "        \n",
    "        # if we have solved at least once and the solution is converged\n",
    "        if gwf.solve(1) and kiter > 1:\n",
    "            # get the sfr dsflow rate for the last reach\n",
    "            addr = [\"DSFLOW\", \"FREYBERG6\", \"SFR_1\"]\n",
    "            wbaddr = gwf.get_var_address(*addr)\n",
    "            sfr_dsflow = gwf.get_value(wbaddr)\n",
    "            print(\"flow stress period,time step {0},{1} converged with {2} iters, final sfr:{3:5.2f}\".\\\n",
    "                  format(stress_period,time_step, kiter,sfr_dsflow[-1]))\n",
    "            break\n",
    "        \n",
    "        \n",
    "        print(\"  ---  iter:{0}, sfr:{1:5.2f} org wel:{2:5.2f} new wel:{3:5.2f} fac:{4:2.1f}\" .\\\n",
    "              format(kiter,sfr_dsflow[-1],org,np.sum(wel_bound),fac))\n",
    "        kiter += 1\n",
    "    try:\n",
    "        gwf.finalize_solve(1)\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    \n",
    "    gwf.finalize_time_step()\n",
    "    # update current sim time\n",
    "    ctime = gwf.get_current_time()\n",
    "gwf.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3477dd5d",
   "metadata": {},
   "source": [
    "Woah - that drastically increased the number of iterations - we are working the solver here with some strong nonlinearity!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ebaab6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "lst_org = flopy.utils.Mf6ListBudget(os.path.join(work_dir,\"freyberg6.lst\"))\n",
    "lst_api = flopy.utils.Mf6ListBudget(os.path.join(api_dir2,\"freyberg6.lst\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e81e2441",
   "metadata": {},
   "outputs": [],
   "source": [
    "oinc,ocum = lst_org.get_dataframes(diff=True)\n",
    "ainc,acum = lst_api.get_dataframes(diff=True)\n",
    "ainc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f22b9e4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axes = plt.subplots(2,1,figsize=(10,10))\n",
    "axes[0].plot(ainc.index,ainc.wel,color=\"r\",label=\"api\")\n",
    "axes[0].plot(oinc.index,oinc.wel,color=\"c\",label=\"org\")\n",
    "axes[1].plot(ainc.index,ainc.sfr,color=\"r\",label=\"api\")\n",
    "axes[1].plot(oinc.index,oinc.sfr,color=\"c\",label=\"org\")\n",
    "axes[0].legend()\n",
    "axes[0].set_title(\"well flux\")\n",
    "axes[1].legend()\n",
    "axes[1].set_title(\"sw-gw flux\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca032c43",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "92b93ee0",
   "metadata": {},
   "source": [
    "We can see the effect of our \"management\" strategy:  substantially more gw is extracted during the wet season while less groundwater is extracted during the dry season compared to the original simulation (that used constant extraction).  And we see that this contributes to less groundwater contribution to surface water during the wet season."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed44308d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}