Initial commit

.github/workflows/gha-tests.yml

@@ -0,0 +1,36 @@
+name: gha-tests
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  test:
+    # Specify the environment to run on Windows
+    runs-on: windows-latest
+
+    strategy:
+      matrix:
+        python-version: ['3.10', '3.11', '3.12']  # Test against multiple Python versions, should correspond to requires-python in pyproject.toml
+
+    steps:
+      # Step 1: Checkout the code
+      - uses: actions/checkout@v2
+
+      # Step 2: Set up Python
+      - name: Set up Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      # Step 3: Install dependencies
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e .
+          pip install pytest
+
+      # Step 4: Run pytest
+      - name: Run tests with pytest
+        run: |
+          pytest

.gitignore

@@ -0,0 +1,6 @@
+**/.ipynb_checkpoints
+**/__pycache__
+/dist
+/adepy.egg-info
+/.vscode
+*.pyc

README.md

@@ -0,0 +1,3 @@
+# adepy
+
+Analytical solutions for solute transport in groundwater with Python 

adepy/__init__.py

@@ -0,0 +1,7 @@
+from ._version import __version__
+
+import adepy.oneD as oneD
+import adepy.twoD as twoD
+import adepy.threeD as threeD
+
+

adepy/_version.py

@@ -0,0 +1 @@
+__version__ = "0.1.0"

adepy/oneD.py

@@ -0,0 +1,122 @@
+from scipy.special import erfc
+from scipy.optimize import brentq
+import numpy as np
+
+def finite1(c0, x, t, v, D, L, lamb=0, nterm=1000):
+
+    x = np.atleast_1d(x)
+    t = np.atleast_1d(t)
+
+    if len(x) > 1 and len(t) > 1:
+        raise ValueError('Either x or t should have length 1')
+
+    Pe = v * L / D
+
+    # find roots
+    def betaf(b):
+        return b * 1 / np.tan(b) + Pe / 2
+
+    intervals = [np.pi * i for i  in range(nterm)]
+    betas = []
+    for i in range(len(intervals) - 1):
+        mi = intervals[i] + 1e-10
+        ma = intervals[i + 1] - 1e-10
+        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)
+
+    if lamb == 0:
+        term0 = 1.0
+    else:
+        u = np.sqrt(v**2 + 4 * lamb * D)
+        term0 = (np.exp((v - u) * x / (2 * D)) + (u - v) / (u + v) * np.exp((v + u) * x / (2 * D) - u * L / D)) /\
+                    (1 + (u - v) / (u + v) * np.exp(-u * L / D))  
+
+    term1 = -2 * np.exp(v * x / (2 * D) - v**2 * t / (4 * D) - lamb * t)
+
+    return c0 * (term0 + term1 * series)
+
+def finite3(c0, x, t, v, D, L, lamb=0, nterm=1000):
+    # https://github.com/BYL4746/columntracer/blob/main/columntracer.py
+
+    x = np.atleast_1d(x)
+    t = np.atleast_1d(t)
+
+    if len(x) > 1 and len(t) > 1:
+        raise ValueError('Either x or t should have length 1')
+
+    Pe = v * L / D
+
+    # find roots
+    def betaf(b):
+        return b * 1 / np.tan(b) - b**2 / Pe + Pe / 4
+
+    intervals = [np.pi * i for i in range(nterm)]
+    betas = []
+    for i in range(len(intervals) - 1):
+        mi = intervals[i] + 1e-10
+        ma = intervals[i + 1] - 1e-10
+        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)
+
+    if lamb == 0:
+        term0 = 1.0
+    else:
+        u = np.sqrt(v**2 + 4 * lamb * D)
+        term0 = (np.exp((v - u) * x / (2 * D)) + (u - v) / (u + v) * np.exp((v + u) * x / (2 * D) - u * L / D)) /\
+                ((u + v) / (2 * v) - (u - v)**2 / (2 * v * (u + v)) * np.exp(-u * L / D)) 
+
+    term1 = -2 * Pe * np.exp(v * x / (2 * D) - v**2 * t / (4 * D) - lamb * t)
+
+    return c0 * (term0 + term1 * series)
+
+def seminf1(c0, x, t, v, D, lamb=0):
+    x = np.atleast_1d(x)
+    t = np.atleast_1d(t)
+
+    u = np.sqrt(v**2 + 4 * lamb * D)
+    term = np.exp(x * (v - u) / (2 * D)) * erfc((x - u * t) / (2 * np.sqrt(D * t))) + \
+            np.exp(x * (v + u) / (2 * D)) * erfc((x + u * t) / (2 * np.sqrt(D * t)))
+
+    return c0 * 0.5 * term
+
+def seminf3(c0, x, t, v, D, lamb=0):
+    x = np.atleast_1d(x)
+    t = np.atleast_1d(t)
+
+    u = np.sqrt(v**2 + 4 * lamb * D)
+    if lamb == 0:
+        term = 0.5 * erfc((x - v * t) / (2 * np.sqrt(D * t))) + np.sqrt(t * v**2 / (np.pi * D)) * np.exp(-(x - v * t)**2 / (4 * D * t)) - \
+            0.5 * (1 + v * x / D + t * v**2 / D) * np.exp(v * x / D) * erfc((x + v * t) / (2 * np.sqrt(D * t)))
+        term0 = 1.0
+    else:
+        term = 2 * np.exp(x * v / D - lamb * t) * erfc((x + v * t) / (2 * np.sqrt(D * t))) +\
+            (u / v - 1) * np.exp(x * (v - u) / (2 * D)) * erfc((x - u * t) / (2 * np.sqrt(D * t))) -\
+            (u / v - 1) * np.exp(x * (v + u) / (2 * D)) * erfc((x + u * t) / (2 * np.sqrt(D * t)))
+        term0 = v**2 / (4 * lamb * D)
+
+    return c0 * term0 * term
+

adepy/threeD.py

@@ -0,0 +1,99 @@
+from scipy.special import erfc
+from scipy.integrate import quad
+import numpy as np
+
+def point3(c0, x, y, z, t, v, n, Dx, Dy, Dz, Q, xc, yc, zc, lamb=0):
+    x = np.atleast_1d(x)
+    y = np.atleast_1d(y)
+    z = np.atleast_1d(z)
+    t = np.atleast_1d(t)
+
+    beta = np.sqrt(v**2 + 4 * Dx * lamb)
+    gamma = np.sqrt((x - xc)**2 + Dx * (y - yc)**2 / Dy + Dx * (z - zc)**2 / Dz)
+
+    a = np.exp(v * (x - xc) / (2 * Dx)) / (8 * n * np.pi * gamma * np.sqrt(Dy * Dz))
+    b = np.exp(gamma * beta / (2 * Dx)) * erfc((gamma + beta * t) / (2 * np.sqrt(Dx * t))) + \
+            np.exp(- gamma * beta / (2 * Dx)) * erfc((gamma - beta * t) / (2 * np.sqrt(Dx * t)))
+
+    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')
+
+    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
+
+        if m == 0:
+            Om = (z2 - z1) / h
+        else:
+            Om = (np.sin(zeta * z2) - np.sin(zeta * z1)) / (m * np.pi)
+
+        for n in range(nterm):
+            eta = n * np.pi / w
+            beta = np.sqrt(v**2 + 4 * Dx * (Dy * eta**2 + Dz * zeta**2 + lamb))
+
+            if m == 0 and n == 0:
+                Lmn = 0.5
+            elif m > 0 and n > 0:
+                Lmn = 2.0
+            else:
+                Lmn = 1.0
+
+            if n == 0:
+                Pn = (y2 - y1) / w
+            else:
+                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)))
+
+            add = Lmn * Om * Pn * np.cos(zeta * z) * 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 patchi(c0, x, y, z, t, v, Dx, Dy, Dz, y1, y2, z1, z2, lamb=0):
+    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
+
+        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
+
+    integrate_vec = np.vectorize(integrate)
+
+    term = integrate_vec(t, x, y, z)
+    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

@@ -0,0 +1,107 @@
+from scipy.special import erfc
+from scipy.integrate import quad
+import numpy as np
+
+def point2(c0, x, y, t, v, Dx, Dy, n, Q, 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)
+    term0 = Q / (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')
+
+    if len(t) > 1:
+        series = np.zeros_like(t, dtype=np.float64)
+    else:
+        series = np.zeros_like(x, dtype=np.float64)
+
+    for n in range(nterm):
+        eta = n * np.pi / w
+        beta = np.sqrt(v**2 + 4 * Dx * (eta**2 * Dy + lamb))
+
+        if n == 0:
+            Ln = 0.5
+            Pn = (y2 - y1) / w
+        else:
+            Ln = 1
+            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)))
+
+        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):
+
+    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
+
+        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)
+    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):
+
+    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)
+    term0 = x * sigma / np.sqrt(8 * np.pi * Dx) * np.exp(v * x / (2 * Dx))
+
+    return c0 * term0 * term