gauss-legendre + numba erfc + tests

adepy/__init__.py

@@ -1,7 +1,7 @@
 from ._version import __version__
 
-import adepy.oneD as oneD
-import adepy.twoD as twoD
-import adepy.threeD as threeD
+from adepy.oneD import finite1, finite3, seminf1, seminf3
+from adepy.twoD import point2, stripf, stripi, gauss
+from adepy.threeD import point3, patchi, patchf
 
 

adepy/_helpers.py

@@ -0,0 +1,42 @@
+from numba import njit, vectorize
+from numba.extending import get_cython_function_address
+import ctypes
+import numpy as np
+from scipy.special import roots_legendre
+from scipy.integrate import quad
+
+addr = get_cython_function_address("scipy.special.cython_special", "__pyx_fuse_1erfc")
+functype = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double)
+erfc_fn = functype(addr)
+
+@vectorize('float64(float64)')
+def _vec_erfc(x):
+    return erfc_fn(x)
+
+@njit
+def _erfc_nb(x):
+    return _vec_erfc(x)
+
+def _integrate(integrand, t, *args, order=100, method='legendre'):
+
+    if method == 'legendre':
+        roots, weights = roots_legendre(order)
+
+        def integrate(t, *args):
+            roots_adj = roots * (t - 0) / 2 + (0 + t) / 2
+            F = integrand(roots_adj, *args).dot(weights) * (t - 0)/2
+            return F
+
+    elif method == 'quadrature':
+
+        def integrate(t, *args):
+            F = quad(integrand, 0, t, args=args.items)
+            return F
+
+    else:
+        raise ValueError('Integration method should be "legendre" or "quadrature"')
+
+    integrate_vec = np.vectorize(integrate)
+    term = integrate_vec(t, *args)
+    return term
+

adepy/oneD.py

@@ -1,6 +1,19 @@
-from scipy.special import erfc
 from scipy.optimize import brentq
 import numpy as np
+from numba import njit, vectorize
+from adepy._helpers import _erfc_nb as erfc
+
+# @njit
+def _bserie_finite1(betas, x, t, Pe, L, D, lamb):
+    bs = 0.0
+    if lamb == 0:
+        for b in betas:
+            bs += b * np.sin(b * x / L) * (b**2 + (Pe / 2)**2) * np.exp(-b**2 * D * t / L**2) /\
+                    ((b**2 + (Pe / 2)**2 + Pe / 2) * (b**2 + (Pe / 2)**2 + lamb / D * L**2))
+    else:
+        for b in betas:
+            bs += b * np.sin(b * x / L) * np.exp(-b**2 * D * t / L**2) / (b**2 + (Pe / 2)**2 + Pe / 2)
+    return bs
 
 def finite1(c0, x, t, v, D, L, lamb=0, nterm=1000):
 
@@ -24,19 +37,8 @@ def betaf(b):
         betas.append(brentq(betaf, mi, ma))
 
     # calculate infinite sum up to nterm terms
-    def bserie(x, t):
-        bs = 0.0
-        if lamb == 0:
-            for b in betas:
-                bs += b * np.sin(b * x / L) * (b**2 + (Pe / 2)**2) * np.exp(-b**2 * D * t / L**2) /\
-                        ((b**2 + (Pe / 2)**2 + Pe / 2) * (b**2 + (Pe / 2)**2 + lamb / D * L**2))
-        else:
-            for b in betas:
-                bs += b * np.sin(b * x / L) * np.exp(-b**2 * D * t / L**2) / (b**2 + (Pe / 2)**2 + Pe / 2)
-        return bs
-
-    bseries_vec = np.vectorize(bserie)
-    series = bseries_vec(x, t)
+    # bseries_vec = np.vectorize(_bserie_finite1)
+    series = _bserie_finite1(betas, x, t, Pe, L, D, lamb)
 
     if lamb == 0:
         term0 = 1.0
@@ -49,6 +51,14 @@ def bserie(x, t):
 
     return c0 * (term0 + term1 * series)
 
+# @njit
+def _bserie_finite3(betas, x, t, Pe, L, D, lamb):
+    bs = 0.0
+    for b in betas:
+        bs += b * (b * np.cos(b * x / L) + (Pe / 2) * np.sin(b * x / L)) / (b**2 + (Pe / 2)**2 + Pe) *\
+                np.exp(-b**2 * D * t / L**2) / (b**2 + (Pe / 2)**2 + lamb * L**2 / D)
+    return bs
+
 def finite3(c0, x, t, v, D, L, lamb=0, nterm=1000):
     # https://github.com/BYL4746/columntracer/blob/main/columntracer.py
 
@@ -72,15 +82,8 @@ def betaf(b):
         betas.append(brentq(betaf, mi, ma))
 
     # calculate infinite sum up to nterm terms
-    def bserie(x, t):
-        bs = 0.0
-        for b in betas:
-            bs += b * (b * np.cos(b * x / L) + (Pe / 2) * np.sin(b * x / L)) / (b**2 + (Pe / 2)**2 + Pe) *\
-                    np.exp(-b**2 * D * t / L**2) / (b**2 + (Pe / 2)**2 + lamb * L**2 / D)
-        return bs
-
-    bseries_vec = np.vectorize(bserie)
-    series = bseries_vec(x, t)
+    # bseries_vec = np.vectorize(_bserie_finite3)
+    series = _bserie_finite3(betas, x, t, Pe, L, D, lamb)
 
     if lamb == 0:
         term0 = 1.0

adepy/threeD.py

@@ -1,6 +1,7 @@
-from scipy.special import erfc
-from scipy.integrate import quad
 import numpy as np
+from numba import njit
+from adepy._helpers import _erfc_nb as erfc
+from adepy._helpers import _integrate as integrate
 
 def point3(c0, x, y, z, t, v, n, Dx, Dy, Dz, Q, xc, yc, zc, lamb=0):
     x = np.atleast_1d(x)
@@ -17,21 +18,14 @@ def point3(c0, x, y, z, t, v, n, Dx, Dy, Dz, Q, xc, yc, zc, lamb=0):
 
     return c0 * Q * a * b
 
-def patchf(c0, x, y, z, t, v, Dx, Dy, Dz, w, h, y1, y2, z1, z2, lamb=0, nterm=50):
-
-    x = np.atleast_1d(x)
-    y = np.atleast_1d(y)
-    z = np.atleast_1d(z)
-    t = np.atleast_1d(t)
-
-    if len(t) > 1 and (len(x) > 1 or len(y) > 1 or len(z) > 1):
-        raise ValueError('If multiple values for t are specified, only one x, y and z value are allowed')
+# @njit
+def _series_patchf(x, y, z, t, v, Dx, Dy, Dz, w, h, y1, y2, z1, z2, lamb, nterm):
 
     if len(t) > 1:
         series = np.zeros_like(t, dtype=np.float64)
     else:
         series = np.zeros_like(x, dtype=np.float64)
-    
+
     for m in range(nterm):
         zeta = m * np.pi / h
 
@@ -64,36 +58,40 @@ def patchf(c0, x, y, z, t, v, Dx, Dy, Dz, w, h, y1, y2, z1, z2, lamb=0, nterm=50
             add = np.where(np.isneginf(add), 0.0, add)
             add = np.where(np.isnan(add), 0.0, add)
             series += add
+
+    return series
 
-    return c0 * series
-
-def patchi(c0, x, y, z, t, v, Dx, Dy, Dz, y1, y2, z1, z2, lamb=0):
+def patchf(c0, x, y, z, t, v, Dx, Dy, Dz, w, h, y1, y2, z1, z2, lamb=0, nterm=50):
+
     x = np.atleast_1d(x)
     y = np.atleast_1d(y)
     z = np.atleast_1d(z)
     t = np.atleast_1d(t)
 
-    def integrate(t, x, y, z):
-        def integrand(tau, x, y, z):
-            zeta = tau
-            ig = 1 / zeta**3 * np.exp(-(v**2 / (4 * Dx) + lamb) * zeta**4 - x**2 / (4 * Dx * zeta**4)) *\
-                     (erfc((y1 - y) / (2 * zeta**2 * np.sqrt(Dy))) - erfc((y2 - y) / (2 * zeta**2 * np.sqrt(Dy)))) *\
-                     (erfc((z1 - z) / (2 * zeta**2 * np.sqrt(Dz))) - erfc((z2 - z) / (2 * zeta**2 * np.sqrt(Dz))))
-
-            # ig = (tau**(-3 / 2)) * np.exp(-(v**2 / (4 * Dx) + lamb) * tau - x**2 / (4 * Dx * tau)) *\
-            #         (erfc((y1 - y) / (2 * np.sqrt(Dy * tau))) - erfc((y2 - y) / (2 * np.sqrt(Dy * tau)))) *\
-            #         (erfc((z1 - z) / (2 * np.sqrt(Dz * tau))) - erfc((z2 - z) / (2 * np.sqrt(Dz * tau))))
-            return ig
+    if len(t) > 1 and (len(x) > 1 or len(y) > 1 or len(z) > 1):
+        raise ValueError('If multiple values for t are specified, only one x, y and z value are allowed')
+
+    series = _series_patchf(x, y, z, t, v, Dx, Dy, Dz, w, h, y1, y2, z1, z2, lamb, nterm)
 
-        F = quad(integrand, 0, t**(1/4), args=(x, y, z), full_output=1)[0]
-        # F = quad(integrand, 0, t, args=(x, y, z), full_output=1)[0]
-        return F
+    return c0 * series
 
-    integrate_vec = np.vectorize(integrate)
+@njit
+def _integrand_patchi(tau, x, y, z, v, Dx, Dy, Dz, y1, y2, z1, z2, lamb):
+    ig = 1 / tau**3 * np.exp(-(v**2 / (4 * Dx) + lamb) * tau**4 - x**2 / (4 * Dx * tau**4)) *\
+                (erfc((y1 - y) / (2 * tau**2 * np.sqrt(Dy))) - erfc((y2 - y) / (2 * tau**2 * np.sqrt(Dy)))) *\
+                (erfc((z1 - z) / (2 * tau**2 * np.sqrt(Dz))) - erfc((z2 - z) / (2 * tau**2 * np.sqrt(Dz))))
 
-    term = integrate_vec(t, x, y, z)
+    return ig
+
+def patchi(c0, x, y, z, t, v, Dx, Dy, Dz, y1, y2, z1, z2, lamb=0, order=100):
+    x = np.atleast_1d(x)
+    y = np.atleast_1d(y)
+    z = np.atleast_1d(z)
+    t = np.atleast_1d(t)
+
+    term = integrate(_integrand_patchi, t**(1/4), x, y, z, v, Dx, Dy, Dz, y1, y2, z1, z2, lamb, order=order, method='legendre')
+
     term0 = x * np.exp(v * x / (2 * Dx)) / (2 * np.sqrt(np.pi * Dx))
-    # term0 = x * np.exp(v * x / (2 * Dx)) / (8 * np.sqrt(np.pi * Dx))
 
     return c0 * term0 * term
 

adepy/twoD.py

@@ -1,39 +1,28 @@
-from scipy.special import erfc
-from scipy.integrate import quad
 import numpy as np
+from numba import njit
+from adepy._helpers import _erfc_nb as erfc
+from adepy._helpers import _integrate as integrate
+
+@njit
+def _integrand_point2(tau, x, y, v, Dx, Dy, xc, yc, lamb):
+    return 1 / tau * np.exp(-(v**2 / (4 * Dx) + lamb) * tau - (x - xc)**2 / (4 * Dx * tau) - (y - yc)**2 / (4 * Dy * tau))
+
+def point2(c0, x, y, t, v, n, Dx, Dy, Qa, xc, yc, lamb=0, order=100):
 
-def point2(c0, x, y, t, v, Dx, Dy, n, Qa, xc, yc, lamb=0):
     x = np.atleast_1d(x)
     y = np.atleast_1d(y)
     t = np.atleast_1d(t)
 
     if len(t) > 1 and (len(x) > 1 or len(y) > 1):
         raise ValueError('If multiple values for t are specified, only one x and y value are allowed')
 
-    def integrate(t, x, y):
-        def integrand(tau, x, y):
-            return 1 / tau * np.exp(-(v**2 / (4 * Dx) + lamb) * tau - (x - xc)**2 / (4 * Dx * tau) - (y - yc)**2 / (4 * Dy * tau))
-
-        F = quad(integrand, 0, t, args=(x, y), full_output=1)[0]
-        return F
-    integrate_vec = np.vectorize(integrate)
-
-    term = integrate_vec(t, x, y)
+    term = integrate(_integrand_point2, t, x, y, v, Dx, Dy, xc, yc, lamb, order=order, method='legendre')
     term0 = Qa / (4 * n * np.pi * np.sqrt(Dx * Dy)) * np.exp(v * (x - xc) / (2 * Dx))
 
     return c0 * term0 * term
 
-def stripf(c0, x, y, t, v, Dx, Dy, y2, y1, w, lamb=0, nterm=100):
-    x = np.atleast_1d(x)
-    y = np.atleast_1d(y)
-    t = np.atleast_1d(t)
-
-    if lamb != 0 and Dy != 0:
-        raise ValueError('Either Dy or lamb should be zero')
-
-    if len(t) > 1 and (len(x) > 1 or len(y) > 1):
-        raise ValueError('If multiple values for t are specified, only one x and y value are allowed')
-
+# @njit
+def _series_stripf(x, y, t, v, Dx, Dy, y2, y1, w, lamb, nterm):
     if len(t) > 1:
         series = np.zeros_like(t, dtype=np.float64)
     else:
@@ -51,57 +40,62 @@ def stripf(c0, x, y, t, v, Dx, Dy, y2, y1, w, lamb=0, nterm=100):
             Pn = (np.sin(eta * y2) - np.sin(eta * y1)) / (n * np.pi)
 
         term = np.exp((x * (v - beta)) / (2 * Dx)) * erfc((x - beta * t) / (2 * np.sqrt(Dx * t))) +\
-                np.exp((x * (v + beta)) / (2 * Dx)) * erfc((x + beta * t) / (2 * np.sqrt(Dx * t)))
+            np.exp((x * (v + beta)) / (2 * Dx)) * erfc((x + beta * t) / (2 * np.sqrt(Dx * t)))
 
         add = Ln * Pn * np.cos(eta * y) * term
 
         add = np.where(np.isneginf(add), 0.0, add)
         add = np.where(np.isnan(add), 0.0, add)
         series += add
 
-    return c0 * series
-
-def stripi(c0, x, y, t, v, Dx, Dy, y2, y1, lamb=0):
+    return series
 
+def stripf(c0, x, y, t, v, Dx, Dy, y2, y1, w, lamb=0, nterm=100):
     x = np.atleast_1d(x)
     y = np.atleast_1d(y)
     t = np.atleast_1d(t)
 
-    def integrate(t, x, y):
-        # error in Wexler, 1992, eq. 91a: denominator of last erfc term should be multiplied by 2. Correct in code below.
-        def integrand(tau, x, y):
-            ig = (tau**(-3 / 2)) * np.exp(-(v**2 / (4 * Dx) + lamb) * tau - x**2 / (4 * Dx * tau)) *\
-                (erfc((y1 - y) / (2 * np.sqrt(Dy * tau))) - erfc((y2 - y) / (2 * np.sqrt(Dy * tau))))
-            return ig
+    if lamb == 0 and Dy == 0:
+        raise ValueError('Either Dy or lamb should be non-zero')
+
+    if len(t) > 1 and (len(x) > 1 or len(y) > 1):
+        raise ValueError('If multiple values for t are specified, only one x and y value are allowed')
+
+    series = _series_stripf(x, y, t, v, Dx, Dy, y2, y1, w, lamb, nterm)
+
+    return c0 * series
 
-        F = quad(integrand, 0, t, args=(x, y), full_output=1)[0]
-        return F
+@njit
+def _integrand_stripi(tau, x, y, v, Dx, Dy, y2, y1, lamb):
+    # error in Wexler, 1992, eq. 91a: denominator of last erfc term should be multiplied by 2. Correct in code below.
+    ig = (tau**(-3 / 2)) * np.exp(-(v**2 / (4 * Dx) + lamb) * tau - x**2 / (4 * Dx * tau)) *\
+        (erfc((y1 - y) / (2 * np.sqrt(Dy * tau))) - erfc((y2 - y) / (2 * np.sqrt(Dy * tau))))
+    return ig
 
-    integrate_vec = np.vectorize(integrate)
+def stripi(c0, x, y, t, v, Dx, Dy, y2, y1, lamb=0, order=100):
 
-    term = integrate_vec(t, x, y)
+    x = np.atleast_1d(x)
+    y = np.atleast_1d(y)
+    t = np.atleast_1d(t)
+
+    term = integrate(_integrand_stripi, t, x, y, v, Dx, Dy, y2, y1, lamb, order=order, method='legendre')
     term0 = x / (4 * np.sqrt(np.pi * Dx)) * np.exp(v * x / (2 * Dx))
 
     return c0 * term0 * term
 
-def gauss(c0, x, y, t, v, Dx, Dy, yc, sigma, lamb=0):
+@njit
+def _integrand_gauss(tau, x, y, v, Dx, Dy, yc, sigma, lamb):
+    num = np.exp(-(v**2 / (4 * Dx) + lamb) * tau - x**2 / (4 * Dx * tau) - ((y - yc)**2) / (4 * (Dy * tau + 0.5 * sigma**2)))
+    denom = (tau**(3 / 2)) * np.sqrt(Dy * tau + 0.5 * sigma**2)
+    return num / denom
+
+def gauss(c0, x, y, t, v, Dx, Dy, yc, sigma, lamb=0, order=100):
 
     x = np.atleast_1d(x)
     y = np.atleast_1d(y)
     t = np.atleast_1d(t)
 
-    def integrate(t, x, y):
-        def integrand(tau, x, y):
-            num = np.exp(-(v**2 / (4 * Dx) + lamb) * tau - x**2 / (4 * Dx * tau) - ((y - yc)**2) / (4 * (Dy * tau + 0.5 * sigma**2)))
-            denom = (tau**(3 / 2)) * np.sqrt(Dy * tau + 0.5 * sigma**2)
-            return num / denom
-
-        F = quad(integrand, 0, t, args=(x, y), full_output=1)[0]
-        return F
-
-    integrate_vec = np.vectorize(integrate)
-
-    term = integrate_vec(t, x, y)
+    term = integrate(_integrand_gauss, t, x, y, v, Dx, Dy, yc, sigma, lamb, order=order, method='legendre')
     term0 = x * sigma / np.sqrt(8 * np.pi * Dx) * np.exp(v * x / (2 * Dx))
 
     return c0 * term0 * term

pyproject.toml

@@ -10,7 +10,7 @@ description = "Analytical solutions for solute transport in groundwater with Pyt
 readme = "README.md"
 license = { file = "LICENSE" }
 requires-python = ">=3.10"
-dependencies = ["numpy", "scipy"]
+dependencies = ["numpy", "scipy", "numba"]
 classifiers = [
         'Programming Language :: Python :: 3 :: Only',
         'Programming Language :: Python :: 3.10',