update doc/readme; Add another example

README.md

@@ -71,6 +71,10 @@ The following functions are provided directly
 * ``grf`` "General Radial Flow" Model from *Barker 1988*
 * ``ext_grf`` the transient extended GRF model
 * ``ext_grf_steady`` the steady extended GRF model
+* ``ext_thiem_tpl`` extended Thiem solution for truncated power laws
+* ``ext_theis_tpl`` extended Theis solution for truncated power laws
+* ``ext_thiem_tpl_3d`` extended Thiem solution in 3D for truncated power laws
+* ``ext_theis_tpl_3d`` extended Theis solution in 3D for truncated power laws
 
 
 ### Laplace Transformation
@@ -84,8 +88,8 @@ inverse laplace-transformation of a given function
 
 ## Requirements
 
-- [NumPy >= 1.13.0](https://www.numpy.org)
-- [SciPy >= 0.19.1](https://www.scipy.org)
+- [NumPy >= 1.14.5](https://www.numpy.org)
+- [SciPy >= 1.1.0](https://www.scipy.org)
 - [pentapy](https://github.com/GeoStat-Framework/pentapy)
 
 
@@ -96,8 +100,8 @@ You can contact us via <info@geostat-framework.org>.
 
 ## License
 
-[GPL][gpl_link] © 2019
+[MIT][mit_link] © 2019
 
-[gpl_link]: https://github.com/GeoStat-Framework/AnaFlow/blob/master/LICENSE
+[mit_link]: https://github.com/GeoStat-Framework/AnaFlow/blob/master/LICENSE
 [ogs5_link]: https://www.opengeosys.org/ogs-5/
 [doc_link]: https://geostat-framework.readthedocs.io/projects/anaflow/en/latest/

anaflow/_version.py

@@ -1,3 +1,3 @@
 # -*- coding: utf-8 -*-
 """Provide a central version"""
-__version__ = "1.0.0.dev0"
+__version__ = "1.0.0.rc1"

docs/source/index.rst

@@ -38,6 +38,10 @@ The following functions are provided directly
 * :any:`grf` "General Radial Flow" Model from *Barker 1988*
 * :py:func:`ext_grf` the transient extended GRF model
 * :any:`ext_grf_steady` the steady extended GRF model
+* :any:`ext_thiem_tpl` extended Thiem solution for truncated power laws
+* :any:`ext_theis_tpl` extended Theis solution for truncated power laws
+* :any:`ext_thiem_tpl_3d` extended Thiem solution in 3D for truncated power laws
+* :any:`ext_theis_tpl_3d` extended Theis solution in 3D for truncated power laws
 
 
 Laplace Transformation
@@ -53,12 +57,12 @@ inverse laplace-transformation of a given function
 Requirements
 ============
 
-- `NumPy >= 1.13.0 <https://www.numpy.org>`_
-- `SciPy >= 0.19.1 <https://www.scipy.org/>`_
+- `NumPy >= 1.14.5 <https://www.numpy.org>`_
+- `SciPy >= 1.1.0 <https://www.scipy.org/>`_
 - `pentapy <https://github.com/GeoStat-Framework/pentapy>`_
 
 
 License
 =======
 
-`GPL <https://github.com/GeoStat-Framework/AnaFlow/blob/master/LICENSE>`_ © 2019
+`MIT <https://github.com/GeoStat-Framework/AnaFlow/blob/master/LICENSE>`_ © 2019

examples/01_call_theis.py

@@ -7,7 +7,7 @@
 time = [10, 100, 1000]
 rad = np.geomspace(0.1, 10)
 
-head = theis(time=time, rad=rad, transmissivity=1e-4, storage=1e-4, rate=-1e-4)
+head = theis(time=time, rad=rad, storage=1e-4, transmissivity=1e-4, rate=-1e-4)
 
 for i, step in enumerate(time):
     plt.plot(rad, head[i], label="Theis(t={})".format(step))

examples/10_convergence_ext_theis_tpl.py

@@ -6,21 +6,19 @@
 time = 1e4                   # time point for steady state
 rad = np.geomspace(0.1, 10)  # radius from the pumping well in [0, 4]
 r_ref = 10.0                 # reference radius
-KG = 1e-3                    # the geometric mean of the transmissivity
-len_scale = 10.0             # correlation length of the log-transmissivity
+KG = 1e-4                    # the geometric mean of the transmissivity
+len_scale = 5.0              # correlation length of the log-transmissivity
 hurst = 0.5                  # hurst coefficient
 var = 0.5                    # variance of the log-transmissivity
 rate = -1e-4                 # pumping rate
 dim = 1.5                    # using a fractional dimension
 
 head1 = ext_thiem_tpl(rad, r_ref, KG, len_scale, hurst, var, dim=dim, rate=rate)
-head2 = ext_theis_tpl(time, rad, 1e-3, KG, len_scale, hurst, var, dim=dim, rate=rate, r_bound=r_ref)
+head2 = ext_theis_tpl(time, rad, 1e-4, KG, len_scale, hurst, var, dim=dim, rate=rate, r_bound=r_ref)
 
 plt.plot(rad, head1, label="Ext Thiem TPL")
 plt.plot(rad, head2, label="Ext Theis TPL (t={})".format(time), linestyle="--")
 
-# plt.plot(rad, head1 / head2, label="ratio")
-
 plt.xlabel("r in [m]")
 plt.ylabel("h in [m]")
 plt.legend()

examples/11_convergence_ext_grf.py

@@ -6,19 +6,18 @@
 time = 1e4                   # time point for steady state
 rad = np.geomspace(0.1, 10)  # radius from the pumping well in [0, 4]
 r_ref = 10.0                 # reference radius
-K = 1e-3                     # the geometric mean of the transmissivity
+K = 1e-4                     # the geometric mean of the transmissivity
 rate = -1e-4                 # pumping rate
 dim = 1.5                    # using a fractional dimension
 
 head1 = ext_grf_steady(rad, r_ref, K, dim=dim, rate=rate)
-head2 = ext_grf(time, rad, [1e-3], [K], [0, r_ref], dim=dim, rate=rate)
-head3 = grf(time, rad, 1e-3, K, dim=dim, rate=rate)
+head2 = ext_grf(time, rad, [1e-4], [K], [0, r_ref], dim=dim, rate=rate)
+head3 = grf(time, rad, 1e-4, K, dim=dim, rate=rate)
+head3 -= head3[-1]  # quasi-steady
 
 plt.plot(rad, head1, label="Ext GRF steady")
 plt.plot(rad, head2, label="Ext GRF (t={})".format(time), linestyle="--")
-plt.plot(rad, head3 - head3[-1], label="GRF quasi-steady (t={})".format(time), linestyle=":")
-
-# plt.plot(rad, head1 / head2, label="ratio")
+plt.plot(rad, head3, label="GRF quasi-steady (t={})".format(time), linestyle=":")
 
 plt.xlabel("r in [m]")
 plt.ylabel("h in [m]")

examples/12_compare_theis_quasi_steady.py

@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+from matplotlib import pyplot as plt
+from anaflow import theis, thiem
+
+
+time = [10, 100, 1000]
+rad = np.geomspace(0.1, 10)
+r_ref = np.full_like(rad, 10.0)
+
+head_ref = theis(time, 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)
+
+for i, step in enumerate(time):
+    label_1 = "Theis quasi steady" if i == 0 else None
+    label_2 = "Theis bounded" if i == 0 else None
+    plt.plot(rad, head1[i], label=label_1, color="C"+str(i), linestyle="--")
+    plt.plot(rad, head2[i], label=label_2, color="C"+str(i))
+
+plt.plot(rad, head3, label="Thiem", color="k", linestyle=":")
+
+plt.xlabel("r in [m]")
+plt.ylabel("h in [m]")
+plt.legend()
+plt.tight_layout()
+plt.show()