Add sol from Neuman2004; Update examples; update DOC

anaflow/__init__.py

@@ -41,6 +41,7 @@
    ext_thiem3D
    ext_theis3D
    ext_theis_tpl
+   neuman2004
 
 Special
 ^^^^^^^
@@ -76,6 +77,7 @@
     ext_thiem3D,
     ext_theis3D,
     ext_theis_tpl,
+    neuman2004,
     grf_disk,
 )
 from anaflow.tools.laplace import get_lap_inv, get_lap
@@ -89,6 +91,8 @@
     "ext_theis2D",
     "ext_thiem3D",
     "ext_theis3D",
+    "ext_theis_tpl",
+    "neuman2004",
     "grf_model",
     "grf_disk",
     "get_lap_inv",

anaflow/flow/__init__.py

@@ -41,6 +41,7 @@
    ext_thiem3D
    ext_theis3D
    ext_theis_tpl
+   neuman2004
 
 Special
 ~~~~~~~
@@ -61,6 +62,7 @@
     ext_thiem3D,
     ext_theis3D,
     ext_theis_tpl,
+    neuman2004,
 )
 from anaflow.flow.special import grf_disk
 
@@ -73,5 +75,6 @@
     "ext_thiem3D",
     "ext_theis3D",
     "ext_theis_tpl",
+    "neuman2004",
     "grf_disk",
 ]

anaflow/flow/heterogeneous.py

@@ -12,6 +12,7 @@
    ext_theis2D
    ext_theis3D
    ext_theis_tpl
+   neuman2004
 """
 # pylint: disable=C0103
 from __future__ import absolute_import, division, print_function
@@ -20,7 +21,12 @@
 from scipy.special import exp1, expi
 
 from anaflow.tools.laplace import get_lap_inv
-from anaflow.tools.special import aniso, specialrange_cut, Shaper
+from anaflow.tools.special import (
+    aniso,
+    specialrange_cut,
+    Shaper,
+    neuman2004_trans,
+)
 from anaflow.tools.mean import annular_hmean
 from anaflow.tools.coarse_graining import (
     T_CG,
@@ -38,6 +44,7 @@
     "ext_theis2D",
     "ext_theis3D",
     "ext_theis_tpl",
+    "neuman2004",
 ]
 
 
@@ -680,25 +687,27 @@ def ext_theis_tpl(
         corralation-length of conductivity-distribution
     hurst: :class:`float`
         Hurst coefficient of the TPL model. Should be in (0, 1).
-    c : :class:`float`, optional
-        Intensity of variation in the TPL model.
-        Is overwritten if sig2 is given.
-        Default: ``1.0``
     sig2: :class:`float` or :any:`None`
         Log-normal-variance of the conductivity-distribution.
         If sig2 is given, c will be calculated accordingly.
         Default: :any:`None`
+    c : :class:`float`, optional
+        Intensity of variation in the TPL model.
+        Is overwritten if sig2 is given.
+        Default: ``1.0``
     e : :class:`float`, optional
         Anisotropy-ratio of the vertical and horizontal corralation-lengths.
         This is only applied in 3 dimensions.
         Default: 1.0
     dim: :class:`float`, optional
         Dimension of space.
         Default: ``2.0``
-    Qw : :class:`float`
+    lat_ext : :class:`float`, optional
+        Thickness of the aquifer (lateral extend).
+        Default: ``1.0``
+    Qw : :class:`float`, optional
         Pumpingrate at the well
-    L : :class:`float`
-        Thickness of the aquifer
+        Default: ``-1e-4``
     struc_grid : :class:`bool`, optional
         If this is set to ``False``, the `rad` and `time` array will be merged
         and interpreted as single, r-t points. In this case they need to have
@@ -846,6 +855,150 @@ def ext_theis_tpl(
     return res
 
 
+###############################################################################
+# transient solution for apparent transmissity from Neuman 1994
+###############################################################################
+
+
+def neuman2004(
+    time,
+    rad,
+    storage,
+    trans_gmean,
+    var,
+    len_scale,
+    rate=-1e-4,
+    struc_grid=True,
+    r_well=0.0,
+    r_bound=np.inf,
+    h_bound=0.0,
+    parts=30,
+    lap_kwargs=None,
+):
+    """
+    The transient solution for the apparent transmissivity from [Neuman2004].
+
+    This solution is build on the apparent transmissivity from Neuman 1994,
+    which represents a mean drawdown in an ensemble of pumping tests in
+    heterogeneous transmissivity fields following an exponential covariance.
+
+    Parameters
+    ----------
+    time : :class:`numpy.ndarray`
+        Array with all time-points where the function should be evaluated
+    rad : :class:`numpy.ndarray`
+        Array with all radii where the function should be evaluated
+    storage : :class:`float`
+        Given storativity of the aquifer
+    trans_gmean : :class:`float`
+        Geometric-mean transmissivity.
+    var : :class:`float`
+        Variance of log-transmissivity.
+    len_scale : :class:`float`
+        Correlation-length of log-transmissivity.
+    rate : :class:`float`, optional
+        Pumpingrate at the well. Default: -1e-4
+    struc_grid : :class:`bool`, optional
+        If this is set to ``False``, the `rad` and `time` array will be merged
+        and interpreted as single, r-t points. In this case they need to have
+        the same shapes. Otherwise a structured r-t grid is created.
+        Default: ``True``
+    r_well : :class:`float`, optional
+        Inner radius of the pumping-well. Default: ``0.0``
+    r_bound : :class:`float`, optional
+        Radius of the outer boundary of the aquifer. Default: ``np.inf``
+    h_bound : :class:`float`, optional
+        Reference head at the outer boundary as well as initial condition.
+        Default: ``0.0``
+    parts : :class:`int`, optional
+        Since the solution is calculated by setting the transmissity to local
+        constant values, one needs to specify the number of partitions of the
+        transmissivity. Default: ``30``
+    lap_kwargs : :class:`dict` or :any:`None` optional
+        Dictionary for :any:`get_lap_inv` containing `method` and
+        `method_dict`. The default is equivalent to
+        ``lap_kwargs = {"method": "stehfest", "method_dict": None}``.
+        Default: :any:`None`
+
+    Returns
+    -------
+    head : :class:`numpy.ndarray`
+        Array with all heads at the given radii and time-points.
+
+    References
+    ----------
+    .. [Neuman2004] Neuman, Shlomo P., Alberto Guadagnini, and Monica Riva.
+       ''Type-curve estimation of statistical heterogeneity.''
+       Water resources research 40.4, 2004
+    """
+    Input = Shaper(time, rad, struc_grid)
+    lap_kwargs = {} if lap_kwargs is None else lap_kwargs
+
+    # check the input
+    if r_well < 0.0:
+        raise ValueError("The wellradius needs to be >= 0")
+    if r_bound <= r_well:
+        raise ValueError(
+            "The upper boundary needs to be greater than the wellradius"
+        )
+    if Input.rad_min < r_well:
+        raise ValueError(
+            "The given radii need to be greater than the wellradius"
+        )
+    if trans_gmean <= 0.0:
+        raise ValueError("The Transmissivity needs to be positiv")
+    if var <= 0.0:
+        raise ValueError("The variance needs to be positiv")
+    if len_scale <= 0.0:
+        raise ValueError("The correlationlength needs to be positiv")
+    if storage <= 0.0:
+        raise ValueError("The Storage needs to be positiv")
+    if parts <= 1:
+        raise ValueError("The numbor of partitions needs to be at least 2")
+
+    # genearte rlast from a given relativ-error to farfield-transmissivity
+    rlast = 2 * len_scale
+    # generate the partition points
+    if rlast > r_well:
+        rpart = specialrange_cut(r_well, r_bound, parts + 1, rlast)
+    else:
+        rpart = np.array([r_well, r_bound])
+    # calculate the harmonic mean transmissivity values within each partition
+    Tpart = annular_hmean(
+        neuman2004_trans,
+        rpart,
+        trans_gmean=trans_gmean,
+        var=var,
+        len_scale=len_scale,
+    )
+
+    # write the paramters in kwargs to use the grf-model
+    kwargs = {
+        "rad": Input.rad,
+        "Qw": rate,
+        "dim": 2,
+        "lat_ext": 1,
+        "rpart": rpart,
+        "Spart": np.full_like(Tpart, storage),
+        "Kpart": Tpart,
+        "Kwell": neuman2004_trans(r_well, trans_gmean, var, len_scale),
+    }
+    kwargs.update(lap_kwargs)
+
+    res = np.zeros((Input.time_no, Input.rad_no))
+    # call the grf-model
+    lap_inv = get_lap_inv(grf_laplace, **kwargs)
+    res[Input.time_gz, :] = lap_inv(Input.time[Input.time_gz])
+    res = Input.reshape(res)
+    if rate > 0:
+        res = np.maximum(res, 0)
+    else:
+        res = np.minimum(res, 0)
+    # add the reference head
+    res += h_bound
+    return res
+
+
 if __name__ == "__main__":
     import doctest
 

anaflow/flow/homogeneous.py

@@ -285,7 +285,7 @@ def grf_model(
     if S <= 0.0:
         raise ValueError("The Storage needs to be positiv")
     if dim <= 0.0 or dim > 3.0:
-        raise ValueError("The dimension needs to be positiv and < 3")
+        raise ValueError("The dimension needs to be positiv and <= 3")
     if lat_ext <= 0.0:
         raise ValueError("The lateral extend needs to be positiv")
 

anaflow/flow/laplace.py

@@ -91,8 +91,8 @@ def grf_laplace(
         Given storativity values for each disk
     Qw : :class:`float`
         Pumpingrate at the well
-    Twell : :class:`float`, optional
-        Transmissivity at the well. Default: ``Tpart[0]``
+    Kwell : :class:`float`, optional
+        Conductivity at the well. Default: ``Kpart[0]``
     cut_off_prec : :class:`float`, optional
         Define a cut-off precision for the calculation to select the disks
         included in the calculation. Default ``1e-20``
@@ -126,8 +126,7 @@ def grf_laplace(
     # initialize the result
     res = np.zeros(s.shape + rad.shape)
     # set the conductivity at the well
-    if Kwell is None:
-        Kwell = Kpart[0]
+    Kwell = Kpart[0] if Kwell is None else float(Kwell)
     # the first sqrt of the diffusivity values
     diff_sr0 = np.sqrt(Spart[0] / Kpart[0])
     # set the general pumping-condtion depending on the well-radius
@@ -219,8 +218,7 @@ def grf_laplace(
             if rpart[-1] < np.inf:
                 Mb[-2, -2] = kv0(Cs[-1] * rpart[-1])
                 Mb[2, -1] = iv0(Cs[-1] * rpart[-1])
-            else:
-                # erase the last row, since X[-1] will be 0
+            else:  # erase the last row, since X[-1] will be 0
                 Mb[0, -1] = 0
                 Mb[1, -1] = 0
 

anaflow/flow/special.py

@@ -24,12 +24,13 @@
 def grf_disk(
     time,
     rad,
-    K_part,
     S_part,
+    K_part,
     R_part,
-    Qw,
+    Qw=-1e-4,
     dim=2,
     lat_ext=1.0,
+    K_well=None,
     struc_grid=True,
     r_well=0.0,
     r_bound=np.inf,
@@ -51,14 +52,16 @@ def grf_disk(
         Array with all time-points where the function should be evaluated
     rad : :class:`numpy.ndarray`
         Array with all radii where the function should be evaluated
-    K_part : :class:`numpy.ndarray`
-        Given conductivity values for each disk
     S_part : :class:`numpy.ndarray`
         Given storativity values for each disk
+    K_part : :class:`numpy.ndarray`
+        Given conductivity values for each disk
     R_part : :class:`numpy.ndarray`
-        Given radii separating the disks
-    Qw : :class:`float`
-        Pumpingrate at the well
+        Given radii separating the disks (excluding r_well and r_bound).
+    Qw : :class:`float`, optional
+        Pumpingrate at the well. Default: -1e-4
+    K_well : :class:`float`, optional
+        Conductivity at the well. Default: ``K_part[0]``
     struc_grid : :class:`bool`, optional
         If this is set to ``False``, the `rad` and `time` array will be merged
         and interpreted as single, r-t points. In this case they need to have
@@ -107,13 +110,14 @@ def grf_disk(
         "Kpart": K_part,
         "dim": dim,
         "lat_ext": lat_ext,
+        "Kwell": K_well,
     }
     kwargs.update(lap_kwargs)
 
     res = np.zeros((Input.time_no, Input.rad_no))
     # call the grf-model
     lap_inv = get_lap_inv(grf_laplace, **kwargs)
-    res[Input.time > 0, :] = lap_inv(Input.time[Input.time > 0])
+    res[Input.time_gz, :] = lap_inv(Input.time[Input.time_gz])
     res = Input.reshape(res)
     if Qw > 0:
         res = np.maximum(res, 0)

anaflow/tools/coarse_graining.py

@@ -83,10 +83,10 @@ def T_CG(rad, TG, sig2, corr, prop=1.6, Twell=None):
     """
     rad = np.squeeze(rad)
 
-    if Twell is not None:
-        chi = np.log(Twell) - np.log(TG)
+    if Twell is None:
+        chi = -sig2 / 2.0  # T_well = T_H
     else:
-        chi = -sig2 / 2.0
+        chi = np.log(Twell) - np.log(TG)
 
     return TG * np.exp(chi / (1.0 + (prop * rad / corr) ** 2))
 

anaflow/tools/mean.py

@@ -391,7 +391,7 @@ def annular_pmean(func, val_arr, p=2.0, ann_dim=2, arg_dict=None, **kwargs):
         func=func,
         val_arr=val_arr,
         f_def=lambda x: np.power(x, p),
-        f_inv=lambda x: np.power(x, 1 / p),
+        f_inv=lambda x: np.power(x, 1.0 / p),
         ann_dim=ann_dim,
         arg_dict=arg_dict,
         **kwargs