New MODFLOW layout

.gitignore

@@ -87,36 +87,12 @@ ipython_config.py
 #   intended to run in multiple environments; otherwise, check them in:
 # .python-version
 
-# pipenv
-#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
-#   However, in case of collaboration, if having platform-specific dependencies or dependencies
-#   having no cross-platform support, pipenv may install dependencies that don't work, or not
-#   install all needed dependencies.
-#Pipfile.lock
-
-# poetry
-#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
-#   This is especially recommended for binary packages to ensure reproducibility, and is more
-#   commonly ignored for libraries.
-#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
-#poetry.lock
-
 # uv
 #   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
 #   This is especially recommended for binary packages to ensure reproducibility, and is more
 #   commonly ignored for libraries.
 uv.lock
 
-# pdm
-#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
-#pdm.lock
-#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
-#   in version control.
-#   https://pdm.fming.dev/latest/usage/project/#working-with-version-control
-.pdm.toml
-.pdm-python
-.pdm-build/
-
 # PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
 __pypackages__/
 
@@ -160,22 +136,11 @@ dmypy.json
 # Cython debug symbols
 cython_debug/
 
-# PyCharm
-#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
-#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
-#  and can be added to the global gitignore or merged into this file.  For a more nuclear
-#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
-#.idea/
-
 /docs/examples/bin
-/docs/examples/pestf_glm
-/docs/examples/pestf_hp
-/docs/examples/pestf_ies
-/docs/examples/pestf_sen
+/docs/examples/pestf_*
 /docs/examples/template
 /docs/examples/model
 /docs/examples/mf_files
-docs/examples/bin/mf6
-docs/examples/data/evap.csv
-docs/examples/data/obs.csv
-docs/examples/data/rain.csv
+/docs/examples/data/*.csv
+*.code-workspace
+docs/examples/mftest/

docs/examples/modflow_uzf_newstyle.ipynb

@@ -0,0 +1,292 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Using a Modflow UZF model as a stressmodel in Pastas\n",
+    "\n",
+    "Please note that the showcased workflow is not how we intend to keep things. Major breaking changes in the next release."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "\n",
+    "import flopy\n",
+    "import pandas as pd\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\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data Groundwater Article"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = ps.load_dataset(\"collenteur_2019\")\n",
+    "head = ds[\"head\"].squeeze()\n",
+    "prec = ds[\"rain\"].squeeze().resample(\"D\").asfreq().fillna(0.0)\n",
+    "evap = ds[\"evap\"].squeeze()\n",
+    "\n",
+    "tmin = pd.Timestamp(\"2003-01-01\")\n",
+    "tmax = pd.Timestamp(\"2018-12-25\")\n",
+    "\n",
+    "ps.plots.series(head, [prec, evap], hist=False);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Standard exponential model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mlexp = ps.Model(head)\n",
+    "mlexp.add_stressmodel(\n",
+    "    ps.RechargeModel(prec=prec, evap=evap, rfunc=ps.Exponential(), name=\"exp\")\n",
+    ")\n",
+    "mlexp.solve(\n",
+    "    tmin=tmin,\n",
+    "    tmax=tmax,\n",
+    ")\n",
+    "mlexp.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## FLEX model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mlflex = ps.Model(head)\n",
+    "rm = ps.RechargeModel(\n",
+    "    prec=prec.multiply(1e3),\n",
+    "    evap=evap.multiply(1e3),\n",
+    "    rfunc=ps.Exponential(),\n",
+    "    recharge=ps.rch.FlexModel(),\n",
+    "    name=\"flex\",\n",
+    ")\n",
+    "mlflex.add_stressmodel(rm)\n",
+    "mlflex.set_parameter(f\"{rm.name}_srmax\", vary=False)\n",
+    "mlflex.solve(tmin=tmin, tmax=tmax)\n",
+    "mlflex.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### UZF model calibrated"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mlmf = ps.Model(head, name=\"uzftest\")\n",
+    "mfmodel = ppmf.ModflowModel(\n",
+    "    model=mlmf,\n",
+    "    exe_name=bindir / \"mf6\",\n",
+    "    sim_ws=Path(\"uzftest\"),\n",
+    "    silent=True,\n",
+    "    solver_kwargs=dict(\n",
+    "        print_option=\"summary\",\n",
+    "        outer_dvclose=3e-2,\n",
+    "        outer_maximum=300,\n",
+    "        under_relaxation=\"dbd\",\n",
+    "        linear_acceleration=\"BICGSTAB\",\n",
+    "        under_relaxation_theta=0.7,\n",
+    "        under_relaxation_kappa=0.08,\n",
+    "        under_relaxation_gamma=0.05,\n",
+    "        under_relaxation_momentum=0.0,\n",
+    "        inner_dvclose=3e-2,\n",
+    "        rcloserecord=\"1000.0 strict\",\n",
+    "        inner_maximum=500,\n",
+    "        relaxation_factor=0.97,\n",
+    "        number_orthogonalizations=2,\n",
+    "        preconditioner_levels=8,\n",
+    "        preconditioner_drop_tolerance=0.001,\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "# shorten the warmup to speed up the modflow calculation somewhat.\n",
+    "mlmf.settings[\"warmup\"] = pd.Timedelta(days=4 * 365)\n",
+    "ghb = ppmf.ModflowGhb()\n",
+    "uzf = ppmf.ModflowUzf(prec, evap, ntrailwaves=35, nwavesets=100)\n",
+    "mfmodel.add_modflow_package([ghb, uzf])\n",
+    "\n",
+    "mlmf.solve(diff_step=1e-4, max_nfev=3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Compare"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ax = mlexp.observations().plot(\n",
+    "    marker=\".\", color=\"k\", linestyle=\"None\", legend=False, figsize=(7.5, 6)\n",
+    ")\n",
+    "mlexp.simulate().plot(ax=ax, label=\"Exponential\")\n",
+    "mlflex.simulate().plot(ax=ax, label=\"FlexModel\")\n",
+    "mlmf.simulate().plot(ax=ax, label=\"UZF\")\n",
+    "ax.legend(ncol=4, loc=\"upper right\")\n",
+    "ax.set_xlim(tmin + pd.Timedelta(days=2 * 365), tmax)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Synthetic data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = ps.load_dataset(\"vonk_2024\")\n",
+    "head = ds[\"heads\"].loc[:, (\"SandyLoam\", \"D85\", \"x00\")].rename(\"usg_head\").iloc[1:]\n",
+    "prec = ds[\"stresses\"].loc[:, \"prec [m/d]\"].rename(\"prec\")\n",
+    "evap = ds[\"stresses\"].loc[:, \"evap [m/d]\"].rename(\"evap\")\n",
+    "\n",
+    "tmin = pd.Timestamp(\"2000-01-01\")\n",
+    "tmax = head.index[-1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### UZF model with parameters we \"know\" from USG"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mlu = ps.Model(head)\n",
+    "mlu.del_constant()\n",
+    "\n",
+    "mf = ppmf.ModflowUzf(\n",
+    "    exe_name=bindir / \"mf6\",\n",
+    "    sim_ws=\"mf_files/mfuzf_usg\",\n",
+    "    head=head,\n",
+    "    ntrailwaves=35,\n",
+    "    nwavesets=100,\n",
+    ")\n",
+    "sm = ppmf.ModflowModel(prec, evap, modflow=mf, name=\"mfsm\")\n",
+    "mlu.add_stressmodel(sm)\n",
+    "\n",
+    "# only parameter we don't know\n",
+    "mlu.set_parameter(f\"{sm.name}_c\", initial=100.0)  # select by hand\n",
+    "# plugin \"true\" parameters from modflow usg model\n",
+    "mlu.set_parameter(f\"{sm.name}_d\", initial=8.5)  # level in canal\n",
+    "mlu.set_parameter(\n",
+    "    f\"{sm.name}_s\", initial=0.271\n",
+    ")  # specific yield from theta_s - theta_fc\n",
+    "mlu.set_parameter(f\"{sm.name}_height\", initial=1.5)  # freeboard\n",
+    "mlu.set_parameter(f\"{sm.name}_vks\", initial=1.0061)  # m/day\n",
+    "mlu.set_parameter(f\"{sm.name}_thtr\", initial=0.065)  # theta_r\n",
+    "mlu.set_parameter(f\"{sm.name}_thts\", initial=0.41, pmax=0.45)  # theta_s\n",
+    "mlu.set_parameter(f\"{sm.name}_eps\", initial=4.954)  # brooks in mfusg\n",
+    "mlu.set_parameter(f\"{sm.name}_extdpfrac\", initial=1.0 / 1.5)  # root depth 1 m\n",
+    "\n",
+    "mlu.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### UZF model calibrated"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with SolveTimer() as st:\n",
+    "    mlu.solve(callback=st.timer, diff_step=0.1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mlu.plot()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pastas-plugins (3.13.5)",
+   "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.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

pastas_plugins/modflow/__init__.py

@@ -2,6 +2,7 @@
 from pastas_plugins.modflow.modflow import (
     ModflowDrn,
     ModflowDrnSto,
+    ModflowGhb,
     ModflowRch,
     ModflowSto,
     ModflowUzf,