Merge pull request #167 from mbakker7/use_ttim_line_elements

Use ttim line elements

Author: dbrakenhoff

docs/02examples/02_wells_rivers_and_recharge.ipynb

@@ -49,9 +49,7 @@
    "outputs": [],
    "source": [
     "# Create basic model elements\n",
-    "ml = tml.ModelMaq(\n",
-    "    kaq=[2, 6, 4], z=[165, 140, 120, 80, 60, 0], c=[2000, 20000], npor=0.3\n",
-    ")\n",
+    "ml = tml.ModelMaq(kaq=[2, 6, 4], z=[165, 140, 120, 80, 60, 0], c=[2000, 20000], npor=0.3)\n",
     "rf = tml.Constant(ml, xr=20000, yr=20000, hr=175, layer=0)\n",
     "p = tml.CircAreaSink(ml, xc=10000, yc=10000, R=15000, N=0.0002, layer=0)\n",
     "w1 = tml.Well(ml, xw=10000, yw=8000, Qw=1000, rw=0.3, layers=0, label=\"well 1\")\n",
@@ -97,10 +95,10 @@
     "]\n",
     "hls4 = [170, np.nan, 166, np.nan, 162, np.nan, np.nan, 156]\n",
     "\n",
-    "ls1 = tml.HeadLineSinkString(ml, xy=xy1, hls=hls1, layers=0)\n",
-    "ls2 = tml.HeadLineSinkString(ml, xy=xy2, hls=hls2, layers=0)\n",
-    "ls3 = tml.HeadLineSinkString(ml, xy=xy3, hls=hls3, layers=0)\n",
-    "ls4 = tml.HeadLineSinkString(ml, xy=xy4, hls=hls4, layers=0)"
+    "ls1 = tml.RiverString(ml, xy=xy1, hls=hls1, layers=0)\n",
+    "ls2 = tml.RiverString(ml, xy=xy2, hls=hls2, layers=0)\n",
+    "ls3 = tml.RiverString(ml, xy=xy3, hls=hls3, layers=0)\n",
+    "ls4 = tml.RiverString(ml, xy=xy4, hls=hls4, layers=0)"
    ]
   },
   {
@@ -185,6 +183,13 @@
     "w3.plotcapzone(hstepmax=50, nt=10, zstart=30, tmax=50 * 365.25, orientation=\"both\")\n",
     "w3.plotcapzone(hstepmax=50, nt=10, zstart=100, tmax=50 * 365.25, orientation=\"both\")"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

docs/02examples/04_horizontal_wells.ipynb

@@ -238,9 +238,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ml.plots.vcontour(\n",
-    "    win=[-400, 400, 0, 0], n=200, levels=20, vinterp=False, figsize=(16, 4)\n",
-    ")"
+    "ml.plots.vcontour(win=[-400, 400, 0, 0], n=200, levels=20, vinterp=False, figsize=(16, 4))"
    ]
   },
   {
@@ -293,9 +291,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ml.plots.vcontour(\n",
-    "    win=[-400, 400, 0, 0], n=200, levels=20, vinterp=False, figsize=(16, 4)\n",
-    ")"
+    "ml.plots.vcontour(win=[-400, 400, 0, 0], n=200, levels=20, vinterp=False, figsize=(16, 4))"
    ]
   },
   {

docs/02examples/building_pit.ipynb

@@ -194,7 +194,7 @@
     "    layers=layers,\n",
     ")\n",
     "\n",
-    "tml.HeadLineSink(\n",
+    "tml.River(\n",
     "    ml,\n",
     "    x1=-length / 2 + offset,\n",
     "    y1=width / 2 - offset,\n",
@@ -203,7 +203,7 @@
     "    hls=h_bem,\n",
     "    layers=np.arange(last_lay_dw + 1),\n",
     ")\n",
-    "tml.HeadLineSink(\n",
+    "tml.River(\n",
     "    ml,\n",
     "    x1=-length / 2 + offset,\n",
     "    y1=0,\n",
@@ -212,7 +212,7 @@
     "    hls=h_bem,\n",
     "    layers=np.arange(last_lay_dw + 1),\n",
     ")\n",
-    "tml.HeadLineSink(\n",
+    "tml.River(\n",
     "    ml,\n",
     "    x1=-length / 2 + offset,\n",
     "    y1=-width / 2 + offset,\n",
@@ -228,7 +228,7 @@
     "Qtot = 0.0\n",
     "\n",
     "for e in ml.elementlist:\n",
-    "    if e.name == \"HeadLineSink\":\n",
+    "    if e.name == \"River\":\n",
     "        Qtot += e.discharge()\n",
     "\n",
     "print(\"\\nDischarge =\", np.round(Qtot.sum(), 2), \"m^3/day\")\n",
@@ -471,7 +471,7 @@
     "wlayers = np.arange(np.sum(-14 <= ml.aq.zaqbot))\n",
     "wlayers = wlayers[1:]\n",
     "\n",
-    "tml.HeadLineSink(\n",
+    "tml.River(\n",
     "    ml,\n",
     "    x1=-length / 2 + offset,\n",
     "    y1=width / 2 - offset,\n",
@@ -480,7 +480,7 @@
     "    hls=h_bem,\n",
     "    layers=wlayers,\n",
     ")\n",
-    "tml.HeadLineSink(\n",
+    "tml.River(\n",
     "    ml,\n",
     "    x1=-length / 2 + offset,\n",
     "    y1=0,\n",
@@ -490,7 +490,7 @@
     "    layers=wlayers,\n",
     "    order=5,\n",
     ")\n",
-    "tml.HeadLineSink(\n",
+    "tml.River(\n",
     "    ml,\n",
     "    x1=-length / 2 + offset,\n",
     "    y1=-width / 2 + offset,\n",
@@ -506,7 +506,7 @@
     "Qtot_ml = 0.0\n",
     "\n",
     "for e in ml.elementlist:\n",
-    "    if e.name == \"HeadLineSink\":\n",
+    "    if e.name == \"River\":\n",
     "        Qtot_ml += e.discharge()\n",
     "\n",
     "print(\"\\nDischarge =\", np.round(Qtot_ml.sum(), 2), \"m^3/day\")\n",

docs/04tests/test_linesink_discharge.ipynb

@@ -46,7 +46,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(\"head at center of line-sink:\", ml1.head(ls1.xc, ls1.yc))\n",
+    "print(\"head at center of line-sink:\", ml1.head(ls1.xc[0], ls1.yc[0]))\n",
     "print(\"discharge of line-sink:\", ls1.discharge())"
    ]
   },
@@ -732,16 +732,28 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ml.plots.vcontour(\n",
-    "    win=[-20, 20, 0, 0], n=100, levels=20, layout=True, horizontal_axis=\"x\"\n",
-    ")\n",
+    "ml.plots.vcontour(win=[-20, 20, 0, 0], n=100, levels=20, layout=True, horizontal_axis=\"x\")\n",
     "plt.xlabel(\"x (m)\");"
    ]
   }
  ],
  "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/04tests/test_normal_flux.ipynb

@@ -223,9 +223,7 @@
     "for i in range(len(xy_out) - 1):\n",
     "    xyi = np.array(xy_out[i : i + 2])\n",
     "    (p2,) = plt.plot(xyi[:, 0], xyi[:, 1], ls=\"dashed\", color=\"r\", label=\"well outside\")\n",
-    "plt.legend(\n",
-    "    [p1, p2], [p1.get_label(), p2.get_label()], loc=(0, 1), frameon=False, ncol=2\n",
-    ")\n",
+    "plt.legend([p1, p2], [p1.get_label(), p2.get_label()], loc=(0, 1), frameon=False, ncol=2)\n",
     "plt.grid(True)"
    ]
   },

docs/04tests/test_polygon_areasink.ipynb

@@ -48,9 +48,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ml.plots.contour(\n",
-    "    win=[-100, 100, -100, 100], ngr=100, layers=[0], levels=20, labels=False\n",
-    ")"
+    "ml.plots.contour(win=[-100, 100, -100, 100], ngr=100, layers=[0], levels=20, labels=False)"
    ]
   },
   {

pyproject.toml

@@ -61,7 +61,7 @@ packages = ["timml", "timml.besselaesnumba"]
 version = { attr = "timml.version.__version__" }
 
 [tool.ruff]
-line-length = 88
+line-length = 90
 extend-include = ["*.ipynb"]
 
 [tool.ruff.lint]

timml/aquifer_parameters.py

@@ -129,9 +129,7 @@ def param_3d(kaq, z, kzoverkh, npor, top="conf", topres=0):
         )
     H = z[:-1] - z[1:]
     assert np.all(H >= 0), "Error: Not all layers thicknesses are non-negative" + str(H)
-    c = 0.5 * H[:-1] / (kzoverkh[:-1] * kaq[:-1]) + 0.5 * H[1:] / (
-        kzoverkh[1:] * kaq[1:]
-    )
+    c = 0.5 * H[:-1] / (kzoverkh[:-1] * kaq[:-1]) + 0.5 * H[1:] / (kzoverkh[1:] * kaq[1:])
     if top == "conf":
         c = np.hstack((1e100, c))
     elif top == "semi":

timml/besselaesnumba/besselaesnumba.py

@@ -479,12 +479,12 @@ def IntegralF(zin, z1in, z2in, Lin, labda, order, Rconv, lstype):
     cln1 = np.log(cd1minz)
     cln2 = np.log(cd2minz)
     for n in range(0, 2 * NTERMS + order + 1):
-        cInt = cInt + (
-            2.0 * calpha[n] * cln2 - 2.0 * calpha[n] / (n + 1) + cbeta[n]
-        ) * (cd2minz) ** (n + 1) / float(n + 1)
-        cInt = cInt - (
-            2.0 * calpha[n] * cln1 - 2.0 * calpha[n] / (n + 1) + cbeta[n]
-        ) * (cd1minz) ** (n + 1) / float(n + 1)
+        cInt = cInt + (2.0 * calpha[n] * cln2 - 2.0 * calpha[n] / (n + 1) + cbeta[n]) * (
+            cd2minz
+        ) ** (n + 1) / float(n + 1)
+        cInt = cInt - (2.0 * calpha[n] * cln1 - 2.0 * calpha[n] / (n + 1) + cbeta[n]) * (
+            cd1minz
+        ) ** (n + 1) / float(n + 1)
     pot = np.real(cInt) * biglab / (2.0 * np.pi)
     return pot
 
@@ -624,12 +624,12 @@ def IntegralG(zin, z1in, z2in, Lin, labda, order, Rconv, lstype):
     cln1 = np.conj(cln1)
     cln2 = np.conj(cln2)
     for n in range(0, 2 * NTERMS + order):
-        g3 = g3 - (calpha[n] * cln2 - calpha[n] / (n + 1)) * (cd2minz) ** (
+        g3 = g3 - (calpha[n] * cln2 - calpha[n] / (n + 1)) * (cd2minz) ** (n + 1) / float(
             n + 1
-        ) / float(n + 1)
-        g3 = g3 + (calpha[n] * cln1 - calpha[n] / (n + 1)) * (cd1minz) ** (
+        )
+        g3 = g3 + (calpha[n] * cln1 - calpha[n] / (n + 1)) * (cd1minz) ** (n + 1) / float(
             n + 1
-        ) / float(n + 1)
+        )
     g3 = biglabin * g3 / (2.0 * np.pi)
     wdis = g1 + g2 + g3
     return wdis

timml/besselaesnumba/besselaesnumba_old.py

@@ -467,12 +467,12 @@ def IntegralF(zin, z1in, z2in, Lin, labda, order, Rconv, lstype):
     cln1 = np.log(cd1minz)
     cln2 = np.log(cd2minz)
     for n in range(0, 2 * NTERMS + order + 1):
-        cInt = cInt + (
-            2.0 * calpha[n] * cln2 - 2.0 * calpha[n] / (n + 1) + cbeta[n]
-        ) * (cd2minz) ** (n + 1) / float(n + 1)
-        cInt = cInt - (
-            2.0 * calpha[n] * cln1 - 2.0 * calpha[n] / (n + 1) + cbeta[n]
-        ) * (cd1minz) ** (n + 1) / float(n + 1)
+        cInt = cInt + (2.0 * calpha[n] * cln2 - 2.0 * calpha[n] / (n + 1) + cbeta[n]) * (
+            cd2minz
+        ) ** (n + 1) / float(n + 1)
+        cInt = cInt - (2.0 * calpha[n] * cln1 - 2.0 * calpha[n] / (n + 1) + cbeta[n]) * (
+            cd1minz
+        ) ** (n + 1) / float(n + 1)
     pot = np.real(cInt) * biglab / (2.0 * np.pi)
     return pot
 
@@ -612,12 +612,12 @@ def IntegralG(zin, z1in, z2in, Lin, labda, order, Rconv, lstype):
     cln1 = np.conj(cln1)
     cln2 = np.conj(cln2)
     for n in range(0, 2 * NTERMS + order):
-        g3 = g3 - (calpha[n] * cln2 - calpha[n] / (n + 1)) * (cd2minz) ** (
+        g3 = g3 - (calpha[n] * cln2 - calpha[n] / (n + 1)) * (cd2minz) ** (n + 1) / float(
             n + 1
-        ) / float(n + 1)
-        g3 = g3 + (calpha[n] * cln1 - calpha[n] / (n + 1)) * (cd1minz) ** (
+        )
+        g3 = g3 + (calpha[n] * cln1 - calpha[n] / (n + 1)) * (cd1minz) ** (n + 1) / float(
             n + 1
-        ) / float(n + 1)
+        )
     g3 = biglabin * g3 / (2.0 * np.pi)
     wdis = g1 + g2 + g3
     return wdis