Update main to v0.2.0

docs/examples/crosscorrelation_seriesJ.ipynb

@@ -25,9 +25,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import pandas as pd\n",
     "import matplotlib as mpl\n",
     "import matplotlib.pyplot as plt\n",
+    "import pandas as pd\n",
     "\n",
     "import pastas_plugins.cross_correlation as ppcc"
    ]
@@ -96,11 +96,17 @@
    ],
    "source": [
     "nlags = 25\n",
-    "ccf_xy = ppcc.ccf(Y, X, alpha=0.05, nlags=nlags) # note the order of the arguments\n",
-    "acf_x = ppcc.ccf(X, X, alpha=0.05, nlags=nlags) # acf is ccf with itself\n",
-    "acf_y = ppcc.ccf(Y, Y, alpha=0.05, nlags=nlags) # acf is ccf with itself\n",
+    "ccf_xy = ppcc.ccf(Y, X, alpha=0.05, nlags=nlags)  # note the order of the arguments\n",
+    "acf_x = ppcc.ccf(X, X, alpha=0.05, nlags=nlags)  # acf is ccf with itself\n",
+    "acf_y = ppcc.ccf(Y, Y, alpha=0.05, nlags=nlags)  # acf is ccf with itself\n",
     "\n",
-    "f, ax = plt.subplots(3, 1, sharex=True, figsize=(6.5, 5.0), constrained_layout=True,)\n",
+    "f, ax = plt.subplots(\n",
+    "    3,\n",
+    "    1,\n",
+    "    sharex=True,\n",
+    "    figsize=(6.5, 5.0),\n",
+    "    constrained_layout=True,\n",
+    ")\n",
     "ppcc.plot_corr(ccf_xy, ax=ax[0])\n",
     "ax[0].set_ylabel(\"CCF\")\n",
     "ax[0].set_title(\"Cross-correlation between input and output\")\n",

docs/examples/responses.ipynb

@@ -16,9 +16,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import numpy as np\n",
     "import pandas as pd\n",
     "import pastas as ps\n",
-    "import numpy as np\n",
     "\n",
     "import pastas_plugins as pp\n",
     "from pastas_plugins import responses"

pastas_plugins/modflow/modflow.py

@@ -1,10 +1,14 @@
+import functools
+import logging
 from typing import List, Protocol
 
 import flopy
 import numpy as np
 from pandas import DataFrame, Series
 from pastas.typing import ArrayLike
 
+logger = logging.getLogger(__name__)
+
 
 class Modflow(Protocol):
     def __init__(self) -> None: ...
@@ -17,15 +21,21 @@ def simulate(self) -> ArrayLike: ...
 
 
 class ModflowRch:
-    def __init__(self, exe_name: str, sim_ws: str):
-        # self.initialhead = initialhead
+    def __init__(
+        self, exe_name: str, sim_ws: str, raise_on_modflow_error: bool = False
+    ):
         self.exe_name = exe_name
         self.sim_ws = sim_ws
         self._name = "mf_rch"
         self._stress = None
         self._simulation = None
         self._gwf = None
-        self._changing_packages = ("IC", "NPF", "STO", "GHB", "RCH")
+        self._changing_packages = (
+            "STO",
+            "GHB",
+            "RCH",
+        )
+        self.raise_on_modflow_error = raise_on_modflow_error
 
     def get_init_parameters(self, name: str) -> DataFrame:
         parameters = DataFrame(columns=["initial", "pmin", "pmax", "vary", "name"])
@@ -40,6 +50,7 @@ def create_model(self) -> None:
             version="mf6",
             exe_name=self.exe_name,
             sim_ws=self.sim_ws,
+            lazy_io=True,
         )
 
         _ = flopy.mf6.ModflowTdis(
@@ -52,25 +63,15 @@ def create_model(self) -> None:
         gwf = flopy.mf6.ModflowGwf(
             sim,
             modelname=self._name,
-            newtonoptions="NEWTON",
         )
 
         _ = flopy.mf6.ModflowIms(
             sim,
-            # print_option="SUMMARY",
             complexity="SIMPLE",
-            # outer_dvclose=1e-9,
-            # outer_maximum=100,
-            # under_relaxation="DBD",
-            # under_relaxation_theta=0.7,
-            # inner_maximum=300,
-            # inner_dvclose=1e-9,
-            # rcloserecord=1e-3,
+            outer_dvclose=1e-2,
+            inner_dvclose=1e-2,
+            rcloserecord=1e-1,
             linear_acceleration="BICGSTAB",
-            # scaling_method="NONE",
-            # reordering_method="NONE",
-            # relaxation_factor=0.97,
-            # filename=f"{name}.ims",
             pname=None,
         )
         # sim.register_ims_package(imsgwf, [self._name])
@@ -83,12 +84,18 @@ def create_model(self) -> None:
             ncol=1,
             delr=1,
             delc=1,
-            top=1000,
-            botm=-1000,
+            top=1.0,
+            botm=0.0,
             idomain=1,
             pname=None,
         )
 
+        _ = flopy.mf6.ModflowGwfnpf(
+            gwf, save_flows=False, icelltype=0, k=1.0, pname="npf"
+        )
+
+        _ = flopy.mf6.ModflowGwfic(gwf, strt=0.0, pname="ic")
+
         _ = flopy.mf6.ModflowGwfoc(
             gwf,
             head_filerecord=f"{gwf.name}.hds",
@@ -103,9 +110,7 @@ def create_model(self) -> None:
     def update_model(self, p: ArrayLike):
         sy, c, f = p[0:3]
 
-        d = 0
-        laytyp = 1
-
+        d = 0.0
         r = self._stress[0] + f * self._stress[1]
 
         # remove existing packages
@@ -114,19 +119,12 @@ def update_model(self, p: ArrayLike):
         ):
             [self._gwf.remove_package(x) for x in self._changing_packages]
 
-        ic = flopy.mf6.ModflowGwfic(self._gwf, strt=d, pname="ic")
-        ic.write()
-
-        npf = flopy.mf6.ModflowGwfnpf(
-            self._gwf, save_flows=False, icelltype=laytyp, pname="npf"
-        )
-        npf.write()
-
+        haq = (self._gwf.dis.top.array - self._gwf.dis.botm.array)[0]
         sto = flopy.mf6.ModflowGwfsto(
             self._gwf,
             save_flows=False,
-            iconvert=laytyp,
-            sy=sy,
+            iconvert=0,
+            ss=sy / haq,
             transient=True,
             pname="sto",
         )
@@ -136,7 +134,7 @@ def update_model(self, p: ArrayLike):
         ghb = flopy.mf6.ModflowGwfghb(
             self._gwf,
             maxbound=1,
-            stress_period_data={0: [[(0, 0, 0), d, 1 / c]]},
+            stress_period_data={0: [[(0, 0, 0), d, 1.0 / c]]},
             pname="ghb",
         )
         ghb.write()
@@ -162,16 +160,27 @@ def update_model(self, p: ArrayLike):
 
         self._gwf.name_file.write()
 
+    @functools.lru_cache(maxsize=5)
+    def _get_head(self, p):
+        self.update_model(p=p)
+        success, _ = self._simulation.run_simulation(silent=True)
+        if success:
+            return self._gwf.output.head().get_ts((0, 0, 0))[:, 1]
+        else:
+            logger.error(
+                "ModflowError: model run failed with parameters: "
+                f"sy={p[0]}, c={p[1]}, f={p[2]}"
+            )
+            if self.raise_on_modflow_error:
+                raise Exception(
+                    "Modflow run failed. Check the LIST file for more information."
+                )
+            else:
+                return np.zeros(self._nper)
+
     def simulate(self, p: ArrayLike, stress: List[Series]) -> ArrayLike:
         if self._simulation is None:
             self._stress = stress
             self._nper = len(self._stress[0])
             self.create_model()
-        self.update_model(p=p)
-        success, _ = self._simulation.run_simulation(silent=True)
-        if success:
-            heads = flopy.utils.HeadFile(
-                self._simulation.sim_path / f"{self._simulation.name}.hds"
-            ).get_ts((0, 0, 0))
-            return heads[:, 1]
-        return np.zeros(self._stress[0].shape)
+        return self._get_head(tuple(p))

pastas_plugins/modflow/stressmodels.py

@@ -64,7 +64,7 @@ def set_init_parameters(self) -> None:
         self.parameters = self.modflow.get_init_parameters(self.name)
 
     def to_dict(self, series: bool = True) -> dict:
-        pass
+        raise NotImplementedError()
 
     def get_stress(
         self,
@@ -101,7 +101,6 @@ def simulate(
             data=h,
             index=stress[0].index,
             name=self.name,
-            fastpath=True,
         )
 
     def _get_block(self, p, dt, tmin, tmax):

pastas_plugins/modflow/version.py

@@ -1 +1 @@
-__version__ = "0.1.0"
+__version__ = "0.2.0"

pastas_plugins/responses/parameters.py

@@ -0,0 +1,130 @@
+from numpy import pi, sqrt
+from scipy.special import k0
+
+
+def kraijenhoff_parameters(
+    S: float, K: float, D: float, x: float, L: float, N: float = 1.0
+) -> tuple[float]:
+    """Get Pastas parameters for Kraijenhoff van de Leur response function for
+    an homogeneous aquifer between two parallel canals.
+
+    Parameters
+    ----------
+    S : float
+        The storativity of the aquifer [-].
+    K : float
+        The saturated hydraulic conductivity [L/T].
+    D : float
+        The saturated thickness of the aquifer [L].
+    x : float
+        The location in the aquifer [L] with x=0 in the center.
+    L : float
+        The aquifer length [L].
+    N : float
+        The recharge flux [L/T].
+
+    Returns
+    -------
+    tuple[float]
+        A tuple containing the response parameters A, a, and b.
+    """
+    b = x / L
+    A = -N * L**2 / (2 * K * D) * (b**2 - (1 / 4))
+    a = S * L**2 / (pi**2 * K * D)
+    return A, a, b
+
+
+def exponential_parameters(S: float, c: float, N: 1.0) -> tuple[float]:
+    """Get Pastas parameters for an Exponential response for a linear
+    reservoir system.
+
+    Parameters
+    ----------
+    S : float
+        The storativity of the aquifer [-].
+    c : float
+        The drainage resistance [T].
+    N : float
+        The recharge flux [L/T].
+
+    Returns
+    -------
+    tuple[float]
+        A tuple containing the response parameters A and a.
+    """
+
+    A = N * c
+    a = c * S
+    return A, a
+
+
+def hantush_parameters(
+    S: float, K: float, D: float, c: float, r: float, Q: float = -1.0
+) -> tuple[float]:
+    """Get Pastas parameters for an Hantush response for a well in a confined aquifer.
+
+    Parameters
+    ----------
+    S : float
+        The storativity of the aquifer [-].
+    K : float
+        The saturated hydraulic conductivity [L/T].
+    D : float
+        The saturated thickness of the aquifer [L].
+    c : float
+        The drainage resistance [T].
+    r : float
+        The distance from the well resistance [L].
+    Q : float
+        The discharge of the well [L].
+
+    Returns
+    -------
+    tuple[float]
+        A tuple containing the response parameters A, a and b.
+    """
+
+    T = K * D
+    lab = sqrt(T * c)
+    b = r**2 / (4.0 * lab**2.0)
+    a = c * S
+    A = Q * k0(r / lab) / (2.0 * pi * T)
+
+    return A, a, b
+
+
+def theis_parameters(
+    S: float, K: float, D: float, x: float, L: float, Q: float = -1.0
+) -> tuple[float]:
+    """Get Pastas parameters for an Theis response for a well between two parallel canals.
+
+    Parameters
+    ----------
+    S : float
+        The storativity of the aquifer [-].
+    K : float
+        The saturated hydraulic conductivity [L/T].
+    D : float
+        The saturated thickness of the aquifer [L].
+    x : float
+        The location in the aquifer [L] with x=0 in the center.
+    L : float
+        The aquifer length [L].
+    Q : float
+        The discharge of the well [L].
+
+    Note
+    ----
+    Works only along the line y=0
+
+    Returns
+    -------
+    tuple[float]
+        A tuple containing the response parameters A, a and b.
+    """
+    T = K * D
+    xw = 0.0  # with xw = 0
+    A = Q / (4.0 * pi * T)
+    a = S * L**2 / (pi**2.0 * T)
+    b = (x - xw) / L
+    return A, a, b