update examples

examples/09_compare_exttheis2d_neuman.py

@@ -23,6 +23,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.xlabel("r in [m]")
 plt.ylabel("h in [m]")
 plt.legend()

examples/12_compare_theis_quasi_steady.py

@@ -6,12 +6,12 @@
 
 time = [10, 100, 1000]
 rad = np.geomspace(0.1, 10)
-r_ref = np.full_like(rad, 10.0)
+r_ref = 10.0
 
-head_ref = theis(time, r_ref, storage=1e-3, transmissivity=1e-4, rate=-1e-4)
+head_ref = theis(time, np.full_like(rad, r_ref), storage=1e-3, transmissivity=1e-4, rate=-1e-4)
 head1 = theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4) - head_ref
-head2 = theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4, r_bound=r_ref[0])
-head3 = thiem(rad, r_ref[0], transmissivity=1e-4, rate=-1e-4)
+head2 = theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4, r_bound=r_ref)
+head3 = thiem(rad, r_ref, transmissivity=1e-4, rate=-1e-4)
 
 for i, step in enumerate(time):
     label_1 = "Theis quasi steady" if i == 0 else None

examples/13_self_defined_transmissivity.py

@@ -0,0 +1,76 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+from matplotlib import pyplot as plt
+import matplotlib.gridspec as gridspec
+from anaflow import ext_grf, ext_grf_steady
+from anaflow.tools import specialrange_cut, annular_hmean
+
+
+def step_f(rad, R_part, K_part):
+    """Step Transmissivity callable."""
+    return np.piecewise(rad, [np.logical_and(r1 <= rad, rad < r2) for r1, r2 in zip(R_part[:-1], R_part[1:])], K_part)
+
+
+def cond(rad, K_far, K_well, len_scale):
+    """Conductivity with linear increase from K_well to K_far."""
+    rad = np.abs(rad, dtype=float)
+    if K_far > K_well:
+        return np.minimum(K_well + rad / len_scale * (K_far - K_well), K_far)
+    return np.maximum(K_well + rad / len_scale * (K_far - K_well), K_far)
+
+
+time_labels = ["10 s", "100 s", "1000 s"]
+time = [10, 100, 1000]
+rad = np.geomspace(0.1, 6)
+S = 1e-4
+K_well = 1e-5
+K_far = 1e-4
+len_scale = 5.0
+rate = -1e-4
+dim = 1.5
+
+cut_off = len_scale
+parts = 30
+r_well = 0.0
+r_bound = 50.0
+
+# calculate a disk-distribution of "trans" by calculating harmonic means
+R_part = specialrange_cut(r_well, r_bound, parts, cut_off)
+K_part = annular_hmean(cond, R_part, ann_dim=dim, K_far=K_far, K_well=K_well, len_scale=len_scale)
+S_part = np.full_like(K_part, S)
+# calculate transient and steady heads
+head1 = ext_grf(time, rad, S_part, K_part, R_part, dim=dim, rate=rate)
+head2 = ext_grf_steady(rad, r_bound, cond, dim=dim, rate=-1e-4, K_far=K_far, K_well=K_well, len_scale=len_scale)
+
+# plotting
+gs = gridspec.GridSpec(2, 1, height_ratios=[1, 3])
+ax1 = plt.subplot(gs[0])
+ax2 = plt.subplot(gs[1], sharex=ax1)
+time_ticks=[]
+for i, step in enumerate(time):
+    label = "Transient" if i == 0 else None
+    ax2.plot(rad, head1[i], label=label, color="C"+str(i))
+    time_ticks.append(head1[i][-1])
+
+ax2.plot(rad, head2, label="Steady", color="k", linestyle=":")
+
+rad_lin = np.linspace(rad[0], rad[-1], 1000)
+ax1.plot(rad_lin, step_f(rad_lin, R_part, K_part), label="step Conductivity")
+ax1.plot(rad_lin, cond(rad_lin, K_far, K_well, len_scale), label="Conductivity")
+ax1.set_yticks([K_well, K_far])
+ax1.set_xticklabels([])
+ax1.set_ylabel(r"$K$ in $[\frac{m}{s}]$")
+ax1.legend()
+
+ax2.set_xlabel("r in [m]")
+ax2.set_ylabel("h in [m]")
+ax2.legend()
+ylim = ax2.get_ylim()
+ax2.set_xlim([0, rad[-1]])
+ax3 = ax2.twinx()
+ax3.set_yticks(time_ticks)
+ax3.set_yticklabels(time_labels)
+ax3.set_ylim(ylim)
+
+plt.tight_layout()
+plt.show()