update examples; integer-model examples

Author: MuellerSeb

examples/09_compare_exttheis2d_neuman.py

@@ -41,7 +41,7 @@
     plt.plot(rad, head2[i], label=label2, color="C" + str(i), linestyle="--")
     time_ticks.append(head1[i][-1])
 
-plt.title("$T_G={}$, $\sigma^2={}$, $\ell={}$, $S={}$".format(TG, var, len_scale, S))
+plt.title(r"$T_G={}$, $\sigma^2={}$, $\ell={}$, $S={}$".format(TG, var, len_scale, S))
 plt.xlabel("r in [m]")
 plt.ylabel("h in [m]")
 plt.legend()

examples/16_integral_model.py

@@ -0,0 +1,56 @@
+r"""
+The integral variogram model
+============================
+"""
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+import anaflow as af
+
+
+def dashes(i=1, max_n=12, width=1):
+    return i * [width, width] + [max_n * 2 * width - 2 * i * width, width]
+
+
+rad = np.linspace(0, 3, 1000)
+rough = [0.01, 0.1, 0.25, 0.5, 1.0, 1.5, 2.0, 10.0]
+# rescale for integral scale of 1
+length = [1 / (nu * np.sqrt(np.pi) / (2 * nu + 2.0)) for nu in rough]
+ln_gau = 2 / np.sqrt(np.pi)
+var = 1
+TG = 1e-4
+TH = TG * np.exp(-var / 2)
+
+for i, (ln, nu) in enumerate(zip(length, rough)):
+    t = 1e4 * af.tools.Int_CG(rad, TG, var, ln, nu)
+    plt.plot(
+        rad,
+        t,
+        label=r"$T^{Int}_{CG}$($\nu$" + f"={nu:.4})",
+        dashes=dashes(i),
+        color="C0" if nu < 1 else "C2",
+    )
+plt.plot(
+    rad,
+    af.tools.T_CG(rad, TG, var, ln_gau) / TG,
+    label=r"$T_{CG}$ - Gaussian",
+    color="C1",
+)
+plt.title(f"$T_G=1$ cm²/s, $\\sigma^2={var}$")
+plt.xlabel(r"r / $\ell$")
+plt.ylabel(r"T / cm²s⁻¹")
+plt.grid()
+ax = plt.gca()
+ax.legend()
+ax.ticklabel_format(axis="y", style="scientific", scilimits=[-3, 0], useMathText=True)
+ax.locator_params(tight=True, nbins=5)
+
+ylim = ax.get_ylim()
+ax2 = ax.twinx()
+ax2.set_ylim(ylim)
+ax2.set_yticks([TH * 1e4, TG * 1e4])
+ax2.set_yticklabels([r"$T_H$", r"$T_G$"])
+
+plt.tight_layout()
+plt.show()

examples/17_compare_nugget.py

@@ -0,0 +1,86 @@
+"""
+Impact of roughness on effective drawdown on anti-persistent fields
+===================================================================
+
+When the field roughness approaches zero, the drawdown approaches
+the homogeneous solution for the geometric mean transmissivity.
+"""
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+import anaflow as af
+
+min_s = 60
+hour_s = min_s * 60
+day_s = hour_s * 24
+
+
+def dashes(i=1, max_n=12, width=1):
+    return i * [width, width] + [max_n * 2 * width - 2 * i * width, width]
+
+
+time_labels = ["10 s", "30 min", "steady\nshape"]
+time = [10, 30 * min_s, 10 * day_s]  # 10s, 30min, 10d
+rad = np.geomspace(0.05, 4)  # radius from the pumping well in [0, 4]
+
+TG = 1e-4  # the geometric mean of the transmissivity
+var = 2.25  # variance of the log-transmissivity
+len_scale = 10  # correlation length of the log-transmissivity
+
+# different roughness values from very rough to very smooth
+rough = [0.01, 0.1, 0.25, 0.5, 1.0]
+S = 1e-4  # storativity
+rate = -1e-4  # pumping rate
+heads = []
+hom_theis = af.theis(time, rad, S, TG, rate=rate)
+ext_theis = af.ext_theis_2d(
+    time, rad, S, TG, var, len_scale / np.sqrt(np.pi) * 2, rate=rate
+)
+
+for roughness in rough:
+    # rescale to constant integral scale
+    ln = len_scale / (roughness * np.sqrt(np.pi) / (2 * roughness + 2.0))
+    # ln = len_scale
+    heads.append(
+        af.ext_theis_int(time, rad, S, TG, var, ln, roughness, rate=rate, parts=50)
+    )
+
+time_ticks = []
+for i, step in enumerate(time):
+    for j, roughness in enumerate(rough):
+        label = (
+            r"Theis$^{Int}_{CG}(\nu=" + f"{roughness:.4})$"
+            if i == len(time) - 1
+            else None
+        )
+        plt.plot(
+            rad,
+            heads[j][i],
+            label=label,
+            color=f"C0",
+            dashes=dashes(j),
+            alpha=(1 + i) / len(time),
+        )
+    time_ticks.append(hom_theis[i][-1])
+    label = "Theis($T_G$)" if i == len(time) - 1 else None
+    plt.plot(rad, hom_theis[i], label=label, color="k")  # , dashes=dashes(1))
+
+
+plt.title(
+    f"$T_G=1$ cm²/s, $\\sigma^2={var}$, $\\ell={len_scale}$ m, $S={S}$, $Q={-rate*1000}$ L/s"
+)
+plt.xlabel("r / m")
+plt.ylabel("h / m")
+plt.grid()
+plt.legend()
+plt.gca().locator_params(tight=True, nbins=5)
+plt.gca().set_ylim([-3.2, 0.1])
+plt.gca().set_xlim([0, rad[-1]])
+ylim = plt.gca().get_ylim()
+ax2 = plt.gca().twinx()
+ax2.set_yticks(time_ticks)
+ax2.set_yticklabels(time_labels)
+ax2.set_ylim(ylim)
+plt.tight_layout()
+plt.show()

examples/18_compare_roughness.py

@@ -0,0 +1,87 @@
+r"""
+Impact of roughness on effective drawdown on persistent fields
+==============================================================
+
+When the field gets smoother, the drawdown approaches
+the effective theis solution for a gaussian correlation structure.
+"""
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+import anaflow as af
+
+min_s = 60
+hour_s = min_s * 60
+day_s = hour_s * 24
+
+
+def dashes(i=1, max_n=12, width=1):
+    return i * 2 * [width] + [max_n * 2 * width - 2 * i * width, width]
+
+
+time_labels = ["10 s", "30 min", "steady\nshape"]
+time = [10, 30 * min_s, 10 * day_s]  # 10s, 30min, 10days
+
+rad = np.geomspace(0.05, 4)  # radius from the pumping well in [0, 4]
+
+TG = 1e-4  # the geometric mean of the transmissivity
+var = 2.25  # 3  # correlation length of the log-transmissivity
+len_scale = 10  # variance of the log-transmissivity
+
+# different roughness values from very rough to very smooth
+rough = [1.0, 1.5, 2.0, 10.0]
+
+S = 1e-4  # storativity
+rate = -1e-4  # pumping rate
+heads = []
+hom_theis = af.theis(time, rad, S, TG, rate=rate)
+ext_theis = af.ext_theis_2d(
+    time, rad, S, TG, var, len_scale / np.sqrt(np.pi) * 2, rate=rate
+)
+
+for roughness in rough:
+    # rescale to constant integral scale
+    ln = len_scale / (roughness * np.sqrt(np.pi) / (2 * roughness + 2.0))
+    heads.append(
+        af.ext_theis_int(time, rad, S, TG, var, ln, roughness, rate=rate, parts=50)
+    )
+
+time_ticks = []
+for i, step in enumerate(time):
+    for j, roughness in enumerate(rough):
+        label = (
+            r"Theis$^{Int}_{CG}(\nu=" + f"{roughness:.4})$"
+            if i == len(time) - 1
+            else None
+        )
+        plt.plot(
+            rad,
+            heads[j][i],
+            label=label,
+            color=f"C2",
+            dashes=dashes(j),
+            alpha=(1 + i) / len(time),
+        )
+    time_ticks.append(heads[-1][i][-1])
+    label = "extended Theis" if i == len(time) - 1 else None
+    plt.plot(rad, ext_theis[i], label=label, color="k")
+
+
+plt.title(
+    f"$T_G=1$ cm²/s, $\\sigma^2={var}$, $\\ell={len_scale}$ m, $S={S}$, $Q={-rate*1000}$ L/s"
+)
+plt.xlabel("r / m")
+plt.ylabel("h / m")
+plt.grid()
+plt.legend()
+plt.gca().locator_params(tight=True, nbins=5)
+plt.gca().set_ylim([-3.2, 0.1])
+plt.gca().set_xlim([0, rad[-1]])
+ylim = plt.gca().get_ylim()
+ax2 = plt.gca().twinx()
+ax2.set_yticks(time_ticks)
+ax2.set_yticklabels(time_labels)
+ax2.set_ylim(ylim)
+plt.tight_layout()
+plt.show()

examples/19_convergence_ext_theis_int.py

@@ -0,0 +1,62 @@
+r"""
+Convergence of the extended Theis solutions for the Integral model
+==================================================================
+
+Here we set an outer boundary to the transient solution, so this condition
+coincides with the references head of the steady solution.
+"""
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+from anaflow import ext_theis_int, ext_thiem_int
+
+time = 86400  # time point for steady state (1 day)
+rad = np.geomspace(0.05, 5) / 5  # radius from the pumping well in [0, 1]
+r_ref = 2  # reference radius
+
+KG = 1e-4  # the geometric mean of the transmissivity
+len_scale = 1  # correlation length of the log-transmissivity
+roughness = 1  # roughness parameter
+var = 1  # variance of the log-transmissivity
+S = 1e-4  # storativity
+rate = -1e-4  # pumping rate
+dim = 2
+
+head1 = ext_thiem_int(rad, r_ref, KG, var, len_scale, roughness, dim=dim, rate=rate)
+head2 = ext_theis_int(
+    time,
+    rad,
+    S,
+    KG,
+    var,
+    len_scale,
+    roughness,
+    dim=dim,
+    rate=rate,
+    r_bound=r_ref,
+)
+
+plt.plot(rad, head1, label=r"Thiem$^{Int}_{CG}$", linewidth=3, color="k")
+plt.plot(
+    rad,
+    head2,
+    label=r"Theis$^{Int}_{CG}$" + f"(t=1 day)",
+    linestyle="--",
+    linewidth=3,
+    color="C1",
+)
+plt.title(
+    f"$T_G=1$ cm²/s, $\\sigma^2={var}$, $\\nu={roughness}$, $S={S}$, $Q=0.1$ L/s, "
+    + r"$r_{ref}=2\ell$"
+)
+
+plt.xlabel(r"r / $\ell$")
+plt.ylabel("h / m")
+plt.grid()
+plt.gca().set_ylim([-1.25, 0])
+plt.gca().set_xlim([0, 1])
+plt.gca().locator_params(tight=True, nbins=5)
+plt.tight_layout()
+plt.legend(loc="lower right")
+plt.show()