Test Modflow API

docs/examples/modflow_api.ipynb

@@ -0,0 +1,348 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "e9915657",
+   "metadata": {},
+   "source": [
+    "# Test MODFLOW API"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b1a23206",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "\n",
+    "import flopy\n",
+    "import pastas as ps\n",
+    "from pastas.timer import SolveTimer\n",
+    "\n",
+    "from pastas_plugins import modflow as ppmf\n",
+    "\n",
+    "bindir = Path(\"bin\")\n",
+    "if not (bindir / \"mf6\").exists():\n",
+    "    bindir.mkdir(parents=True, exist_ok=True)\n",
+    "    flopy.utils.get_modflow(bindir, repo=\"modflow6\")\n",
+    "\n",
+    "dll = bindir / \"libmf6.so\"\n",
+    "\n",
+    "ps.set_use_cache(True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2e2f3afc",
+   "metadata": {},
+   "source": [
+    "## Load Data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "23509e17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = ps.load_dataset(\"collenteur_2019\")\n",
+    "head = ds[\"head\"].squeeze().dropna()\n",
+    "prec = ds[\"rain\"].squeeze().dropna().resample(\"D\").asfreq().fillna(0.0)\n",
+    "evap = ds[\"evap\"].squeeze().dropna()\n",
+    "\n",
+    "prec = prec.loc[\"2002-11-01\":]\n",
+    "evap = evap.loc[\"2002-11-01\":]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ca39ddb5",
+   "metadata": {},
+   "source": [
+    "## MODFLOW API"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f68bf1eb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml1 = ps.Model(head, name=\"mftest\")\n",
+    "mfml = ppmf.ModflowModelApi(\n",
+    "    model=ml1,\n",
+    "    dll=dll,\n",
+    "    sim_ws=Path(\"mftest\"),\n",
+    "    silent=False,\n",
+    ")\n",
+    "rch = ppmf.ModflowRch(prec, evap)\n",
+    "ghb = ppmf.ModflowGhb()\n",
+    "mfml.add_modflow_package([rch, ghb])\n",
+    "ml1.add_stressmodel(mfml)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6072a7c6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml1.parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cdab95cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml1.get_parameters(\"mfapi\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1d2839e3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = tuple(ml1.parameters.initial.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2e9366bd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %%time\n",
+    "sim = ml1.simulate(p)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "419a05f4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with SolveTimer() as t:\n",
+    "    ml1.solve(diff_step=1e-4, callback=t.timer)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "43ed7e38",
+   "metadata": {},
+   "source": [
+    "## Classic MF6\n",
+    "\n",
+    "Solve with classic MF6 implementation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e9dabe61",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml2 = ps.Model(head, name=\"mftest\")\n",
+    "mfml2 = ppmf.ModflowModel(\n",
+    "    model=ml2,\n",
+    "    exe_name=bindir / \"mf6\",\n",
+    "    sim_ws=Path(\"mftest\"),\n",
+    "    silent=True,\n",
+    ")\n",
+    "rch = ppmf.ModflowRch(prec, evap)\n",
+    "ghb = ppmf.ModflowGhb()\n",
+    "mfml2.add_modflow_package([rch, ghb])\n",
+    "ml2.add_stressmodel(mfml2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "479682d2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = tuple(ml2.parameters.initial.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bedd27fd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "sim2 = ml2.simulate(p)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d5a13cda",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with SolveTimer() as t:\n",
+    "    ml2.solve(diff_step=1e-4, callback=t.timer)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1d067c0e",
+   "metadata": {},
+   "source": [
+    "## Pastas"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8e58c97e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml3 = ps.Model(head, name=\"pastas\")\n",
+    "rm = ps.RechargeModel(prec, evap, rfunc=ps.Exponential(), name=\"rch\")\n",
+    "ml3.add_stressmodel(rm)\n",
+    "ml3.solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "55eaee2c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ax = ml3.plot(figsize=(8, 3))\n",
+    "sim1 = ml1.simulate()\n",
+    "ax.plot(sim1.index, sim1, label=f\"API ({ml1.stats.rsq():.2%})\")\n",
+    "sim2 = ml2.simulate()\n",
+    "ax.plot(sim2.index, sim2, label=f\"MF6 ({ml2.stats.rsq():.2%})\")\n",
+    "ax.legend(loc=(0, 1), frameon=False, ncol=4)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6ac86f71",
+   "metadata": {},
+   "source": [
+    "## MODFLOW API objects (for testing)"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "33663dab",
+   "metadata": {
+    "vscode": {
+     "languageId": "raw"
+    }
+   },
+   "source": [
+    "import modflowapi\n",
+    "\n",
+    "mf6 = modflowapi.ModflowApi(\"bin/libmf6.so\", working_directory=\"mftest\")\n",
+    "mf6.initialize()\n",
+    "\n",
+    "# get API simulation and model objects\n",
+    "sim = modflowapi.extensions.ApiSimulation.load(mf6)\n",
+    "model = sim.get_model(\"mfapi\")\n",
+    "\n",
+    "# prepare first time step\n",
+    "dt = mf6.get_time_step()\n",
+    "mf6.prepare_time_step(dt)\n",
+    "\n",
+    "# now you can access packages and other solution information\n",
+    "# note the kernel will crash if you do not properly make modflow finish its\n",
+    "# solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5c8b27bf",
+   "metadata": {},
+   "source": [
+    "## MODFLOW with callback"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "52f09582",
+   "metadata": {
+    "vscode": {
+     "languageId": "raw"
+    }
+   },
+   "source": [
+    "def callback_function(sim, callback_step):\n",
+    "    d = mfml.parameters.loc[\"constant_d\", \"initial\"]\n",
+    "    C = mfml.parameters.loc[\"GHB_C\", \"initial\"]\n",
+    "    f = mfml.parameters.loc[\"RCH_f\", \"initial\"]      \n",
+    "    S = mfml.parameters.loc[\"STO_S\", \"initial\"]\n",
+    "    ml = sim.mfsm\n",
+    "\n",
+    "    if callback_step == Callbacks.stress_period_start:\n",
+    "        if sim.kper == 0:\n",
+    "            print(\"Starting stress period 1\")\n",
+    "\n",
+    "            # # GHB\n",
+    "            df = ml.ghb.stress_period_data.dataframe\n",
+    "            df.loc[0, \"bhead\"] = d\n",
+    "            df.loc[0, \"cond\"] = C\n",
+    "            ml.ghb.stress_period_data.dataframe = df\n",
+    "\n",
+    "            # STO\n",
+    "            sto = ml.sto\n",
+    "            sto.set_array(\"ss\", np.array([[[S / 11.0]]]))\n",
+    "\n",
+    "            # DIS\n",
+    "            dis = ml.dis\n",
+    "            dis.set_array(\"bot\", np.array([[[d - 10]]]))\n",
+    "            dis.set_array(\"top\", np.array([[[d + 1.0]]]))\n",
+    "\n",
+    "            # IC\n",
+    "            ic = ml.ic\n",
+    "            ic.set_array(\"strt\", np.array([[[d]]]))\n",
+    "\n",
+    "            # # RCH\n",
+    "            # df = ml.rch.stress_period_data.dataframe\n",
+    "            # df.loc[0, \"recharge\"] = prec.values + f * evap.values\n",
+    "            # ml.rch.stress_period_data.dataframe = df"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pastas-plugins",
+   "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.13.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

pastas_plugins/modflow/__init__.py

@@ -7,5 +7,5 @@
     ModflowSto,
     ModflowUzf,
 )
-from pastas_plugins.modflow.stressmodels import ModflowModel
+from pastas_plugins.modflow.stressmodels import ModflowModel, ModflowModelApi
 from pastas_plugins.modflow.version import __version__