Fix numpy >= 2.4 issues

docs/steady/02examples/circular_area_sink.ipynb

@@ -52,7 +52,7 @@
    "source": [
     "qx = np.zeros_like(x)\n",
     "for i in range(len(x)):\n",
-    "    qx[i], qy = ml.disvec(x[i], 1e-6)\n",
+    "    qx[i : i + 1], qy = ml.disvec(x[i], 1e-6)\n",
     "plt.plot(x, qx)\n",
     "qxb = N * np.pi * R**2 / (2 * np.pi * R)\n",
     "plt.axhline(qxb, color=\"r\", ls=\"--\")\n",

docs/steady/02examples/connected_wells.ipynb

@@ -28,9 +28,19 @@
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "\n",
+    "import timflow\n",
     "import timflow.steady as tfs"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "timflow.show_versions()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -568,8 +578,7 @@
     "    levels=levels,\n",
     "    decimals=2,\n",
     "    layers=[ilay],\n",
-    "    newfig=False,\n",
-    "    linestyles=\"dashed\",\n",
+    "    linestyles=\"solid\",\n",
     "    color=\"C1\",\n",
     ")\n",
     "plt.plot(ws.xcp, ws.ycp, \"kx\");"
@@ -718,7 +727,7 @@
     "    decimals=2,\n",
     "    layers=[ilay],\n",
     "    newfig=False,\n",
-    "    linestyles=\"dashed\",\n",
+    "    linestyles=\"solid\",\n",
     "    color=\"C1\",\n",
     ")\n",
     "plt.plot(ws.xcp, ws.ycp, \"kx\");"

docs/transient/00userguide/tutorials/tutorial0_start_a_model.ipynb

@@ -35,7 +35,6 @@
     "ls = tft.RiverString(ml, xy=[(-25, -30), (-15, 10), (10, 20)], tsandh=\"fixed\", layers=0)\n",
     "w = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n",
     "\n",
-    "ml.initialize()\n",
     "ml.solve(silent=True)"
    ]
   },
@@ -130,6 +129,13 @@
    "source": [
     "ax = ml.plots.head_along_line(-30, 10, 0, 0, t=10, layers=[0, 1, 2])"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

docs/transient/02examples/line_sink_well_sol.ipynb

@@ -105,7 +105,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.figure(figsize=(16, 4))\n",
+    "plt.figure(figsize=(12, 4))\n",
     "for ilay in [0, 1, 2]:\n",
     "    ax = plt.subplot(1, 3, ilay + 1)\n",
     "    ml.plots.contour(\n",
@@ -118,7 +118,8 @@
     "        labels=True,\n",
     "        decimals=2,\n",
     "        ax=ax,\n",
-    "    )"
+    "    )\n",
+    "    ax.set_aspect(\"equal\", adjustable=\"box\")"
    ]
   },
   {

docs/transient/02examples/meandering_river.ipynb

@@ -77,8 +77,22 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "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,

docs/transient/05benchmarks/well_benchmarks.ipynb

@@ -181,7 +181,7 @@
    "outputs": [],
    "source": [
     "def volume(r, t=1):\n",
-    "    return -2 * np.pi * r * ml.head(r, 0, t) * ml.aq.Scoefaq[0]\n",
+    "    return (-2 * np.pi * r * ml.head(r, 0, t) * ml.aq.Scoefaq[0]).squeeze()\n",
     "\n",
     "\n",
     "quad(volume, 1e-5, np.inf)"

pyproject.toml

@@ -97,7 +97,7 @@ addopts = [
     "--cov-report=term",
     "--cov-report=lcov:./coverage/lcov.info",
     "--ignore=docs/steady/02examples/vertical_anisotropy.ipynb",
-    "--ignore=docs/steady/04benchmarks/benchmarking_besselaes.ipynb",
+    "--ignore=docs/steady/04benchmarks/besselnumba_timing.ipynb",
 ]
 markers = [
     "notebooks: run notebooks",

timflow/steady/constant.py

@@ -170,7 +170,7 @@ def equation(self):
                 # else:
                 #    mat[0, ieq:ieq+e. nunknowns] += -1
                 else:
-                    rhs[0] -= e.potentiallayers(
+                    rhs[0:1] -= e.potentiallayers(
                         self.xc[icp], self.yc[icp], self.layers
                     ).sum(0)
         return mat, rhs

timflow/steady/linedoublet.py

@@ -408,17 +408,15 @@ def potinf(self, x, y, aq=None):
         rv = np.zeros((self.Nld, self.ldlist[0].nparam, aq.naq))
         for i in range(self.Nld):
             rv[i] = self.ldlist[i].potinf(x, y, aq)
-        rv.shape = (self.nparam, aq.naq)
-        return rv
+        return rv.reshape((self.nparam, aq.naq))
 
     def disvecinf(self, x, y, aq=None):
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
         rv = np.zeros((2, self.Nld, self.ldlist[0].nparam, aq.naq))
         for i in range(self.Nld):
             rv[:, i] = self.ldlist[i].disvecinf(x, y, aq)
-        rv.shape = (2, self.nparam, aq.naq)
-        return rv
+        return rv.reshape((2, self.nparam, aq.naq))
 
     def plot(self, ax=None, layer=None):
         if ax is None:

timflow/steady/linesink.py

@@ -597,7 +597,7 @@ def initialize(self):
         self.resfac = (
             np.tile(self.res / self.whfac, self.ncp) * self.strengthinf
         )  # this is resfach !
-        self.resfac.shape = (self.ncp, self.nlayers, self.nunknowns)
+        self.resfac = self.resfac.reshape((self.ncp, self.nlayers, self.nunknowns))
         if len(self.hls) == 1:
             self.hc = self.hls * np.ones(self.nparam)
         elif len(self.hls) == self.ncp:  # head specified at control points
@@ -707,7 +707,7 @@ def equation(self):
                     mat[0, ieq : ieq + e.nunknowns] = self.dischargeinf()
                     break
                 ieq += e.nunknowns
-        rhs[0] = self.Qls
+        rhs[0:1] = self.Qls
         return mat, rhs
 
     def setparams(self, sol):
@@ -740,14 +740,16 @@ def __init__(
         self.xy = np.atleast_2d(xy).astype("d")
         if closed:
             self.xy = np.vstack((self.xy, self.xy[0]))
-        self.x, self.y = self.xy[:, 0], self.xy[:, 1]
         self.order = order  # same for all segments in string
         self.dely = dely  # same for all segments in string
         self.lslist = []
         if self.xy.shape[1] == 2:
             self.nls = len(self.xy) - 1
+            self.x, self.y = self.xy[:, 0], self.xy[:, 1]
         elif self.xy.shape[1] == 4:
             self.nls = len(self.xy)
+            self.x = None
+            self.y = None
         if self.layers.ndim == 1:
             if len(self.layers) == self.nls:
                 self.layers = self.layers[:, np.newaxis]
@@ -801,8 +803,7 @@ def potinf(self, x, y, aq=None):
         if aq in self.aq:
             for i, ls in enumerate(self.lslist):
                 rv[i] = ls.potinf(x, y, aq)
-        rv.shape = (self.nparam, aq.naq)
-        return rv
+        return rv.reshape((self.nparam, aq.naq))
 
     def disvecinf(self, x, y, aq=None):
         if aq is None:
@@ -811,8 +812,7 @@ def disvecinf(self, x, y, aq=None):
         if aq in self.aq:
             for i, ls in enumerate(self.lslist):
                 rv[:, i] = ls.disvecinf(x, y, aq)
-        rv.shape = (2, self.nparam, aq.naq)
-        return rv
+        return rv.reshape((2, self.nparam, aq.naq))
 
     def dischargeinf(self):
         rv = np.zeros((self.nls, self.lslist[0].nparam))
@@ -829,7 +829,7 @@ def discharge_per_linesink(self):
             array of shape (nlay, nlinesinks)
         """
         Qls = self.parameters[:, 0] * self.dischargeinf()
-        Qls.shape = (self.nls, self.nlayers, self.order + 1)
+        Qls = Qls.reshape((self.nls, self.nlayers, self.order + 1))
         Qls = Qls.sum(axis=2)
         rv = np.zeros((self.model.aq.naq, self.nls))
         for i, q in enumerate(Qls):
@@ -840,7 +840,7 @@ def discharge(self):
         """Discharge of the element in each layer."""
         rv = np.zeros(self.aq[0].naq)
         Qls = self.parameters[:, 0] * self.dischargeinf()
-        Qls.shape = (self.nls, self.nlayers, self.order + 1)
+        Qls = Qls.reshape((self.nls, self.nlayers, self.order + 1))
         Qls = np.sum(Qls, 2)
         for i, q in enumerate(Qls):
             rv[self.layers[i]] += q
@@ -947,7 +947,6 @@ def __init__(
         self.res = res
         self.wh = wh
         self.model.add_element(self)
-        # TO DO: TEST FOR DIFFERENT AQUIFERS AND LAYERS
 
     def initialize(self):
         if len(self.hls) == 1:  # one value
@@ -978,6 +977,8 @@ def initialize(self):
                 y1, y2 = self.y[i : i + 2]
             elif self.xy.shape[1] == 4:
                 x1, y1, x2, y2 = self.xy[i]
+            else:
+                raise ValueError("xy should have shape (n, 2) or (n, 4)")
             self.lslist.append(
                 River(
                     self.model,
@@ -1018,7 +1019,7 @@ def equation(self):
         # include resistance by computing position of coefficients in matrix
         # and subtracting resistance terms
         iself = self.model.elementlist.index(self)
-        jcol = np.sum(e.nunknowns for e in self.model.elementlist[:iself])
+        jcol = np.sum([e.nunknowns for e in self.model.elementlist[:iself]]).astype(int)
         irow = 0
         for ls in self.lslist:
             for icp in range(ls.ncp):
@@ -1114,7 +1115,7 @@ def equation(self):
                     mat[0, ieq : ieq + self.nunknowns] = self.dischargeinf()
                     break
                 ieq += e.nunknowns
-        rhs[0] = self.Qls
+        rhs[0:1] = self.Qls
         return mat, rhs
 
 

timflow/steady/well.py

@@ -82,7 +82,7 @@ def initialize(self):
         self.parameters[:, 0] = self.Qw
         self.resfac = self.res / (2 * np.pi * self.rw * self.aq.Haq[self.layers])
         self.resfac = self.resfac * np.identity(self.nlayers)
-        self.resfac.shape = (
+        self.resfac = self.resfac.reshape(
             self.ncp,
             self.nlayers,
             self.nlayers,  # changed to nlayers from nunknowns
@@ -451,7 +451,7 @@ def equation(self):
                     mat[0, ieq : ieq + e.nunknowns] = 1.0
                     break
                 ieq += e.nunknowns
-        rhs[0] = self.Qw
+        rhs[0:1] = self.Qw
         return mat, rhs
 
     def setparams(self, sol):
@@ -618,7 +618,7 @@ def equation(self):
             rhs[i] -= rhs[0]
         # first equation is head at control point equals hcp
         mat[0] = 0.0
-        rhs[0] = self.hcp
+        rhs[0:1] = self.hcp
         aq = self.model.aq.find_aquifer_data(self.xcp, self.ycp)
         ieq = 0
         for e in self.model.elementlist:
@@ -628,7 +628,9 @@ def equation(self):
                 )
                 ieq += e.nunknowns
             else:
-                rhs[0] -= e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp]
+                rhs[0:1] -= (
+                    e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp]
+                )
         return mat, rhs
 
     def setparams(self, sol):
@@ -727,7 +729,7 @@ def potinf(self, x, y, aq=None):
             if r < self.rw:
                 r = self.rw  # If at well, set to at radius
             if aq.ilap:
-                pot[0] = np.log(r / self.rw) / (2 * np.pi)
+                pot[0:1] = np.log(r / self.rw) / (2 * np.pi)
                 pot[1:] = -k0(r / aq.lab[1:]) / (2 * np.pi) / k0(self.rw / aq.lab[1:])
             else:
                 pot[:] = -k0(r / aq.lab) / (2 * np.pi) / k0(self.rw / aq.lab)
@@ -751,8 +753,8 @@ def disvecinf(self, x, y, aq=None):
                 xminxw = self.rw
                 yminyw = 0.0
             if aq.ilap:
-                qx[0] = -1 / (2 * np.pi) * xminxw / rsq
-                qy[0] = -1 / (2 * np.pi) * yminyw / rsq
+                qx[0:1] = -1 / (2 * np.pi) * xminxw / rsq
+                qy[0:1] = -1 / (2 * np.pi) * yminyw / rsq
                 kone = k1(r / aq.lab[1:]) / k0(self.rw / aq.lab[1:])
                 qx[1:] = -kone * xminxw / (r * aq.lab[1:]) / (2 * np.pi)
                 qy[1:] = -kone * yminyw / (r * aq.lab[1:]) / (2 * np.pi)
@@ -846,7 +848,7 @@ def potinf(self, x, y, aq=None):
             for w in self.wlist:
                 rv[j : j + w.nparam] = w.potinf(x, y, aq)
                 j += w.nparam
-        # rv.shape = (self.nparam, aq.naq)
+        # rv = rv.reshape((self.nparam, aq.naq))
         return rv
 
     def disvecinf(self, x, y, aq=None):
@@ -907,7 +909,7 @@ def equation(self):
         # include resistance by finding position of element in coefficient matrix
         # and subtracting resfac (resistance factor).
         iself = self.model.elementlist.index(self)
-        jcol = np.sum(e.nunknowns for e in self.model.elementlist[:iself])
+        jcol = np.sum([e.nunknowns for e in self.model.elementlist[:iself]])
         irow = 0
         for w in self.wlist:
             mat[irow : irow + w.nlayers, jcol : jcol + w.nunknowns] -= w.resfac[
@@ -1004,7 +1006,7 @@ def equation(self):
                     mat[0, ieq : ieq + self.nunknowns] = 1.0
                     break
                 ieq += e.nunknowns
-        rhs[0] = self.Qw
+        rhs[0:1] = self.Qw
         return mat, rhs
 
     def setparams(self, sol):
@@ -1166,7 +1168,7 @@ def equation(self):
 
         # first equation is head at control point equals hcp
         mat[0] = 0
-        rhs[0] = self.hcp
+        rhs[0:1] = self.hcp
         aq = self.model.aq.find_aquifer_data(self.xcp, self.ycp)
         ieq = 0
         for e in self.model.elementlist:
@@ -1184,7 +1186,9 @@ def equation(self):
                     )
                     ieq += e.nunknowns
             else:
-                rhs[0] -= e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp]
+                rhs[0:1] -= (
+                    e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp]
+                )
         return mat, rhs
 
     def setparams(self, sol):

timflow/transient/aquifer.py

@@ -135,10 +135,12 @@ def initialize(self):
             self.coef[:, :, i] = np.linalg.solve(v, b).T
         self.lab = 1.0 / np.sqrt(self.eigval)
         self.lab2 = self.lab.copy()
-        self.lab2.shape = (self.naq, self.model.nint, self.model.npint)
+        self.lab2 = self.lab2.reshape((self.naq, self.model.nint, self.model.npint))
         self.lababs = np.abs(self.lab2[:, :, 0])  # used to check distances
         self.eigvec2 = self.eigvec.copy()
-        self.eigvec2.shape = (self.naq, self.naq, self.model.nint, self.model.npint)
+        self.eigvec2 = self.eigvec2.reshape(
+            (self.naq, self.naq, self.model.nint, self.model.npint)
+        )
 
     def compute_lab_eigvec(self, p, returnA=False, B=None):
         sqrtpSc = np.sqrt(p * self.Scoefll * self.c)