From e95b3f5c03af85fe3dd57358a29d1bfc1e3d3b68 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 18 Feb 2021 15:15:46 +0100 Subject: [PATCH 1/9] Switch to GH actions for CI --- .github/workflows/main.yml | 89 ++++++++++++++++++++++++++++++++++++++ .travis.yml | 81 ---------------------------------- 2 files changed, 89 insertions(+), 81 deletions(-) create mode 100644 .github/workflows/main.yml delete mode 100755 .travis.yml diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml new file mode 100644 index 0000000..1f097b3 --- /dev/null +++ b/.github/workflows/main.yml @@ -0,0 +1,89 @@ +name: Continuous Integration + +on: + push: + branches: + - "master" + - "develop" + tags: + - "*" + pull_request: + branches: + - "develop" + # Allows you to run this workflow manually from the Actions tab + workflow_dispatch: + +env: + # needed by coveralls + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + +jobs: + build_sdist: + name: sdist on ${{ matrix.os }} with py ${{ matrix.python-version }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, windows-latest, macos-latest] + python-version: [3.5, 3.6, 3.7, 3.8, 3.9] + + steps: + - uses: actions/checkout@v2 + with: + fetch-depth: '0' + + - name: Set up Python ${{ matrix.python-version }} + uses: actions\setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -U wheel + pip install -r requirements_setup.txt + pip install -r requirements.txt + pip install -r requirements_test.txt + pip install coveralls>=3.0.0 + + - name: Build sdist + run: | + python setup.py sdist --formats=gztar bdist_wheel -d dist + + - name: Run tests + run: | + python -m pytest --cov anaflow --cov-report term-missing -v tests/ + python -m coveralls --service=github + + - uses: actions/upload-artifact@v2 + if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9' + with: + path: dist + + upload_to_pypi: + needs: [build_sdist] + runs-on: ubuntu-latest + + steps: + - uses: actions/download-artifact@v2 + with: + name: artifact + path: dist + + - name: Publish to Test PyPI + # only if working on develop + if: github.ref == 'refs/heads/master' || github.ref == 'refs/heads/develop' + uses: pypa/gh-action-pypi-publish@master + with: + user: __token__ + password: ${{ secrets.test_pypi_password }} + repository_url: https://test.pypi.org/legacy/ + skip_existing: true + + - name: Publish to PyPI + # only if tagged + if: startsWith(github.ref, 'refs/tags') + uses: pypa/gh-action-pypi-publish@master + with: + user: __token__ + password: ${{ secrets.pypi_password }} diff --git a/.travis.yml b/.travis.yml deleted file mode 100755 index e34b9ec..0000000 --- a/.travis.yml +++ /dev/null @@ -1,81 +0,0 @@ -language: python -python: 3.8 - -# setuptools-scm needs all tags in order to obtain a proper version -git: - depth: false - -env: - global: - # Note: TWINE_PASSWORD is set in Travis settings - - TWINE_USERNAME=geostatframework - - CIBW_BUILD="cp35-* cp36-* cp37-* cp38-*" - # update setuptools to latest version - - CIBW_BEFORE_BUILD="pip install -U setuptools" - # testing with cibuildwheel - - CIBW_TEST_REQUIRES=pytest - - CIBW_TEST_COMMAND="pytest -v {project}/tests" - -notifications: - email: - recipients: - - info@geostat-framework.org - -before_install: - - | - if [[ "$TRAVIS_OS_NAME" = windows ]]; then - choco install python --version 3.8.0 - export PATH="/c/Python38:/c/Python38/Scripts:$PATH" - # make sure it's on PATH as 'python3' - ln -s /c/Python38/python.exe /c/Python38/python3.exe - fi - -install: - - python3 -m pip install cibuildwheel==1.3.0 - -script: - - python3 -m cibuildwheel --output-dir tmp_dist - -stages: - - test - - coverage - - name: deploy - if: (NOT type IN (pull_request)) AND (repo = GeoStat-Framework/AnaFlow) - -jobs: - include: - - stage: test - name: Test on Linux - services: docker - - stage: test - name: Test on MacOS - os: osx - language: generic - - stage: test - name: Test on Windows - os: windows - language: shell - - - stage: coverage - name: Coverage on Linux - services: docker - install: python3 -m pip install .[test] coveralls - script: - - python3 -m pytest --cov anaflow --cov-report term-missing -v tests/ - - python3 -m coveralls - - # Test Deploy source distribution - - stage: deploy - name: Test Deploy - install: python3 -m pip install -U setuptools wheel twine - script: python3 setup.py sdist --formats=gztar bdist_wheel - after_success: - - python3 -m twine upload --verbose --skip-existing --repository-url https://test.pypi.org/legacy/ dist/* - - # Deploy source distribution - - stage: deploy - name: Deploy to PyPI - if: tag IS present - install: python3 -m pip install -U setuptools wheel twine - script: python3 setup.py sdist --formats=gztar bdist_wheel - after_success: python3 -m twine upload --verbose --skip-existing dist/* From a7c56cdadf90d652f80bc1e1d38d3687d0365a63 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 18 Feb 2021 15:15:59 +0100 Subject: [PATCH 2/9] Update README AND LICENSE --- LICENSE | 2 +- README.md | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/LICENSE b/LICENSE index 006289e..51e42f7 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ The MIT License (MIT) -Copyright (c) 2019 - 2020 Sebastian Mueller +Copyright (c) 2019 - 2021 Sebastian Mueller Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index 2479aed..5038433 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,8 @@ [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1135723.svg)](https://doi.org/10.5281/zenodo.1135723) [![PyPI version](https://badge.fury.io/py/anaflow.svg)](https://badge.fury.io/py/anaflow) -[![Build Status](https://travis-ci.com/GeoStat-Framework/AnaFlow.svg?branch=master)](https://travis-ci.com/GeoStat-Framework/AnaFlow) -[![Documentation Status](https://readthedocs.org/projects/docs/badge/?version=stable)](https://anaflow.readthedocs.io/en/stable/) +[![Build Status](https://github.com/GeoStat-Framework/AnaFlow/workflows/Continuous%20Integration/badge.svg?branch=develop)](https://github.com/GeoStat-Framework/AnaFlow/actions) +[![Documentation Status](https://readthedocs.org/projects/docs/badge/?version=latest)](https://anaflow.readthedocs.io/en/latest/) [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/ambv/black)

@@ -25,7 +25,7 @@ You can install the latest version with the following command: ## Documentation for AnaFlow -You can find the documentation under [geostat-framework.readthedocs.io][doc_link]. +You can find the documentation under [https://anaflow.readthedocs.io][doc_link]. ### Example @@ -100,7 +100,7 @@ You can contact us via . ## License -[MIT][mit_link] © 2019 - 2020 +[MIT][mit_link] © 2019 - 2021 [mit_link]: https://github.com/GeoStat-Framework/AnaFlow/blob/master/LICENSE [doc_link]: https://anaflow.readthedocs.io From f26a9cce07acfcad48dd40fad98fbaf353c1e341 Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Tue, 15 Jun 2021 19:03:36 +0200 Subject: [PATCH 3/9] Made the Neuman sloution from Moenche et al. 2001 --- anaflow/flow/Neuman.py | 162 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 anaflow/flow/Neuman.py diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py new file mode 100644 index 0000000..0589571 --- /dev/null +++ b/anaflow/flow/Neuman.py @@ -0,0 +1,162 @@ +# -*- coding: utf-8 -*- +""" +Anaflow subpackage providing the Neuman equation for homogeneous aquifer. + +.. currentmodule:: anaflow.flow.Neuman + +The following functions are provided + +.. autosummary:: + +""" +# pylint: disable=C0103 +import numpy as np +from scipy.special import k0, k1 +from scipy.integrate import quad as integ + +import matplotlib.pyplot as plt + +from anaflow.tools.laplace import get_lap_inv +from anaflow.flow.laplace import grf_laplace +from anaflow.tools.special import Shaper, sph_surf + +__all__ = [] + + +def neuman( + s, + rad, + storage, + transmissivity, + sat_thickness=52.0, + rate=-1e-4, + d=0.6, + l=12.6, + rw=0.1, + kz=0.14, + kr=0.23, + ss=4.27e-5, + n_numbers=52, +): + """ + The Neuman solution. + + + Parameters + ---------- + s : :class:`numpy.ndarray` + Array with all time-points where the function should be evaluated in the Laplace space. + rad : :class:`numpy.ndarray` + Array with all radii where the function should be evaluated. + storage : :class:`float` + Storage of the aquifer. + transmissivity : :class:`float` + Geometric-mean transmissivity. + sat_thickness : :class:`float`, optional + Saturated thickness of the aquifer. + rate : :class:`float`, optional + Pumpingrate at the well. Default: -1e-4 + d : :class:`float`, optional + Vertical distance from initial water table to top of + pumped well screen + l : :class:`float`, optional + Vertical distance from initial water table to bottom + of pumped well screen + rw : :class:`float`, optional + Outside radius of the pumped well screen + kz : :class:`float`, optional + Hydraulic conductivity in the vertical direction + kr : :class:`float`, optional + Hydraulic conductivity in the horizontal direction + Ss : :class:`float`, optional + specific storage + """ + z = float(sat_thickness - l) + kd = kz / kr + s = np.squeeze(s).reshape(-1) + rad = np.squeeze(rad).reshape(-1) + res = np.zeros(s.shape + rad.shape) + + for si, se in enumerate(s): + for ri, re in enumerate(rad): + if re < np.inf: + for n in range(n_numbers): + epsilon_n = np.pi / 2 + n * np.pi + betaw = kd * rw ** 2 + qn = ((epsilon_n ** 2) * betaw + se) ** 0.5 + rd = re / rw + E = k0(qn * rd) * np.cos(epsilon_n * (z / sat_thickness)) / ( + ( + np.sin(epsilon_n * (1 - (d / sat_thickness))) + - np.sin(epsilon_n * (1 - (l / sat_thickness))) + ) + / (qn * k1(qn) * (epsilon_n + 0.5 * np.sin(2 * epsilon_n))) + ) + A = k0(qn) * ( + ( + np.sin(epsilon_n * (1 - (d / sat_thickness))) + - np.sin(epsilon_n * (1 - (l / sat_thickness))) + ) + ** 2 + / ( + epsilon_n + * qn + * k1(qn) + * (epsilon_n + 0.5 * np.sin(2 * epsilon_n)) + ) + ) + E = 2 * E + A = (2 / ((l - d) / sat_thickness)) * A + sw = storage + wd = np.pi * rw ** 2 / (2 * np.pi * rw ** 2 * ss * (l - d)) + res[si, ri] = (2 * E) / (se * ((l - d) / sat_thickness) * (1 + wd * se * (A + sw))) + res *= rate / (2 * np.pi * transmissivity) + return res + + +def neuman_from_laplace( + time, + rad, + storage, + transmissivity, + rate=-1e-4, + h_bound=0.0, + struc_grid=True, + lap_kwargs=None, +): + """Neuman solution form laplace solution.""" + Input = Shaper(time, rad, struc_grid) + lap_kwargs = {} if lap_kwargs is None else lap_kwargs + + if not transmissivity > 0.0: + raise ValueError("The Transmissivity needs to be positive.") + if not storage > 0.0: + raise ValueError("The Storage needs to be positive.") + kwargs = { + "rad": rad, + "storage": storage, + "transmissivity": transmissivity, + "rate": rate, + } + kwargs.update(lap_kwargs) + res = np.zeros((Input.time_no, Input.rad_no)) + lap_inv = get_lap_inv(neuman, **kwargs) + # call the laplace inverse function (only at time points > 0) + res[Input.time_gz, :] = lap_inv(Input.time[Input.time_gz]) + # reshaping results + res = Input.reshape(res) + # add the reference head + res += h_bound + return res + + +time = [10, 600, 36000] # 10s, 10min, 10h +rad = np.geomspace(0.1, 10) +# default parameters +T = 1e-4 +S = 1e-4 +Q = -1e-4 + +neu = neuman_from_laplace(time, rad, storage=S, transmissivity=T) +for i in range(len(time)): + plt.plot(rad, neu[i, :], color="C0", linestyle=":", linewidth=4) From b618afd37e0ee155eac01c64689efbd255c3b6b5 Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Mon, 28 Jun 2021 19:03:01 +0200 Subject: [PATCH 4/9] Made the Neuman sloution from Moenche et al. 1996 --- anaflow/flow/Neuman.py | 166 ++++++++++++++++++++--------------------- 1 file changed, 79 insertions(+), 87 deletions(-) diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py index 0589571..5fe450d 100644 --- a/anaflow/flow/Neuman.py +++ b/anaflow/flow/Neuman.py @@ -7,36 +7,31 @@ The following functions are provided .. autosummary:: - + neuman_unconfined_laplace + neuman_unconfined """ # pylint: disable=C0103 import numpy as np -from scipy.special import k0, k1 -from scipy.integrate import quad as integ - -import matplotlib.pyplot as plt - +from scipy.special import k0 from anaflow.tools.laplace import get_lap_inv -from anaflow.flow.laplace import grf_laplace -from anaflow.tools.special import Shaper, sph_surf +from anaflow.tools.special import Shaper +from scipy.optimize import fsolve __all__ = [] -def neuman( - s, - rad, - storage, - transmissivity, - sat_thickness=52.0, - rate=-1e-4, - d=0.6, - l=12.6, - rw=0.1, - kz=0.14, - kr=0.23, - ss=4.27e-5, - n_numbers=52, +def neuman_unconfined_laplace( + s, + rad, + storage, + transmissivity, + rate, + sat_thickness=52.0, + screen_size=11.88, + well_depth=12.6, + kd=0.61, + specific_yield=0.26, + n_numbers=10, ): """ The Neuman solution. @@ -45,7 +40,8 @@ def neuman( Parameters ---------- s : :class:`numpy.ndarray` - Array with all time-points where the function should be evaluated in the Laplace space. + Array with all "Laplace-space-points" where the function should + be evaluated in the Laplace space. rad : :class:`numpy.ndarray` Array with all radii where the function should be evaluated. storage : :class:`float` @@ -56,23 +52,19 @@ def neuman( Saturated thickness of the aquifer. rate : :class:`float`, optional Pumpingrate at the well. Default: -1e-4 - d : :class:`float`, optional - Vertical distance from initial water table to top of - pumped well screen - l : :class:`float`, optional + screen_size : :class:`float`, optional + Vertical length of the observation screen + well_depth : :class:`float`, optional Vertical distance from initial water table to bottom of pumped well screen - rw : :class:`float`, optional - Outside radius of the pumped well screen - kz : :class:`float`, optional - Hydraulic conductivity in the vertical direction - kr : :class:`float`, optional - Hydraulic conductivity in the horizontal direction - Ss : :class:`float`, optional - specific storage + kd : :class:`float`, optional + Dimensionless parameter for the conductivity. + Kz/Kr : vertical conductivity divided by horizontal conductivity + specific_yield: :class:`float`, optional + specific yield """ - z = float(sat_thickness - l) - kd = kz / kr + z = sat_thickness - well_depth + d = well_depth - screen_size s = np.squeeze(s).reshape(-1) rad = np.squeeze(rad).reshape(-1) res = np.zeros(s.shape + rad.shape) @@ -80,49 +72,61 @@ def neuman( for si, se in enumerate(s): for ri, re in enumerate(rad): if re < np.inf: - for n in range(n_numbers): - epsilon_n = np.pi / 2 + n * np.pi - betaw = kd * rw ** 2 - qn = ((epsilon_n ** 2) * betaw + se) ** 0.5 - rd = re / rw - E = k0(qn * rd) * np.cos(epsilon_n * (z / sat_thickness)) / ( - ( - np.sin(epsilon_n * (1 - (d / sat_thickness))) - - np.sin(epsilon_n * (1 - (l / sat_thickness))) - ) - / (qn * k1(qn) * (epsilon_n + 0.5 * np.sin(2 * epsilon_n))) - ) - A = k0(qn) * ( - ( - np.sin(epsilon_n * (1 - (d / sat_thickness))) - - np.sin(epsilon_n * (1 - (l / sat_thickness))) + rd = re / sat_thickness + betaw = kd * (rd ** 2) + righthand = se / ( + ((storage / specific_yield) * betaw) + se / 1e9 + ) + if righthand > 10: + epsilon_n = np.arange( + np.pi / 2, (n_numbers + 3) * np.pi, np.pi + ) + else: + f = lambda eps: eps * np.tan(eps) - righthand + eps_0 = fsolve( + f, np.arange(np.pi / 2, (n_numbers + 3) * np.pi, np.pi) + ) + epsilon_n = eps_0[0:n_numbers] + print(epsilon_n) + for n in range(min(n_numbers, len(epsilon_n))): + xn = (betaw * (epsilon_n[n] ** 2) + se) ** 0.5 + if epsilon_n[n] == 0: + continue + else: + res[si, ri] = ( + 2 + * k0(xn) + * ( + np.sin( + epsilon_n[n] * (1 - (d / sat_thickness)) + ) + - np.sin( + epsilon_n[n] + * (1 - (well_depth / sat_thickness)) + ) ) - ** 2 - / ( - epsilon_n - * qn - * k1(qn) - * (epsilon_n + 0.5 * np.sin(2 * epsilon_n)) + * np.cos(epsilon_n[n] * (z / sat_thickness)) + ) / ( + se + * ((well_depth - d) / sat_thickness) + * ( + 0.5 * epsilon_n[n] + + 0.25 * np.sin(2 * epsilon_n[n]) ) - ) - E = 2 * E - A = (2 / ((l - d) / sat_thickness)) * A - sw = storage - wd = np.pi * rw ** 2 / (2 * np.pi * rw ** 2 * ss * (l - d)) - res[si, ri] = (2 * E) / (se * ((l - d) / sat_thickness) * (1 + wd * se * (A + sw))) + ) res *= rate / (2 * np.pi * transmissivity) return res -def neuman_from_laplace( - time, - rad, - storage, - transmissivity, - rate=-1e-4, - h_bound=0.0, - struc_grid=True, - lap_kwargs=None, +def neuman_unconfined( + time, + rad, + storage, + transmissivity, + rate, + h_bound=0.0, + struc_grid=True, + lap_kwargs=None, ): """Neuman solution form laplace solution.""" Input = Shaper(time, rad, struc_grid) @@ -140,7 +144,7 @@ def neuman_from_laplace( } kwargs.update(lap_kwargs) res = np.zeros((Input.time_no, Input.rad_no)) - lap_inv = get_lap_inv(neuman, **kwargs) + lap_inv = get_lap_inv(neuman_unconfined_laplace, **kwargs) # call the laplace inverse function (only at time points > 0) res[Input.time_gz, :] = lap_inv(Input.time[Input.time_gz]) # reshaping results @@ -148,15 +152,3 @@ def neuman_from_laplace( # add the reference head res += h_bound return res - - -time = [10, 600, 36000] # 10s, 10min, 10h -rad = np.geomspace(0.1, 10) -# default parameters -T = 1e-4 -S = 1e-4 -Q = -1e-4 - -neu = neuman_from_laplace(time, rad, storage=S, transmissivity=T) -for i in range(len(time)): - plt.plot(rad, neu[i, :], color="C0", linestyle=":", linewidth=4) From 6201b0a8e98a658ce921fab6df190db58220674e Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Wed, 7 Jul 2021 16:44:38 +0200 Subject: [PATCH 5/9] Made the Neuman sloution from Moenche et al. 1996 --- anaflow/flow/Neuman.py | 72 ++++++++++++++++++------------------------ 1 file changed, 31 insertions(+), 41 deletions(-) diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py index 5fe450d..b3cb937 100644 --- a/anaflow/flow/Neuman.py +++ b/anaflow/flow/Neuman.py @@ -16,10 +16,32 @@ from anaflow.tools.laplace import get_lap_inv from anaflow.tools.special import Shaper from scipy.optimize import fsolve +from scipy import optimize __all__ = [] +def eps_func(eps,value=0): + return eps * np.tan(eps) - value + + +def interval_of_nth_root(n, v): + """Correct interval.""" + a, b = n * np.pi, (n + 0.5) * np.pi - 0.001 + if v < 0: # just for completeness (shift to negative branch) + a += 0.5 * np.pi + b += 0.5 * np.pi + return a, b + + +def nth_root(n, v): + """Get nth root of x*tan(x) = v.""" + a, b = interval_of_nth_root(n, v) + if abs(v) > 100: + return (n+0.5) * np.pi + return optimize.bisect(eps_func, a, b, (v,)) + + def neuman_unconfined_laplace( s, rad, @@ -31,7 +53,7 @@ def neuman_unconfined_laplace( well_depth=12.6, kd=0.61, specific_yield=0.26, - n_numbers=10, + n_numbers=300, ): """ The Neuman solution. @@ -74,46 +96,14 @@ def neuman_unconfined_laplace( if re < np.inf: rd = re / sat_thickness betaw = kd * (rd ** 2) - righthand = se / ( - ((storage / specific_yield) * betaw) + se / 1e9 - ) - if righthand > 10: - epsilon_n = np.arange( - np.pi / 2, (n_numbers + 3) * np.pi, np.pi - ) - else: - f = lambda eps: eps * np.tan(eps) - righthand - eps_0 = fsolve( - f, np.arange(np.pi / 2, (n_numbers + 3) * np.pi, np.pi) - ) - epsilon_n = eps_0[0:n_numbers] - print(epsilon_n) - for n in range(min(n_numbers, len(epsilon_n))): - xn = (betaw * (epsilon_n[n] ** 2) + se) ** 0.5 - if epsilon_n[n] == 0: - continue - else: - res[si, ri] = ( - 2 - * k0(xn) - * ( - np.sin( - epsilon_n[n] * (1 - (d / sat_thickness)) - ) - - np.sin( - epsilon_n[n] - * (1 - (well_depth / sat_thickness)) - ) - ) - * np.cos(epsilon_n[n] * (z / sat_thickness)) - ) / ( - se - * ((well_depth - d) / sat_thickness) - * ( - 0.5 * epsilon_n[n] - + 0.25 * np.sin(2 * epsilon_n[n]) - ) - ) + righthand = se / (((storage/specific_yield) * betaw) + se / 1e9) + roots = [nth_root(n, righthand) for n in range(n_numbers)] + for i_n in range(len(roots)): + epsilon_n = roots[i_n] + xn = (betaw * (epsilon_n ** 2) + se) ** 0.5 + res[si, ri] = (2 * k0(xn) * (np.sin(epsilon_n * (1 - (d / sat_thickness))) - np.sin(epsilon_n * + (1 - (well_depth / sat_thickness)))) * np.cos(epsilon_n * (z / sat_thickness))) / \ + (se * ((well_depth - d) / sat_thickness) * (0.5 * epsilon_n + 0.25 * np.sin(2 * epsilon_n))) res *= rate / (2 * np.pi * transmissivity) return res From 6ec28c53de162a45fea4f11191d77f7e82abcf40 Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Fri, 9 Jul 2021 15:03:45 +0200 Subject: [PATCH 6/9] Made the destinction between fully and partially penetrating wells --- anaflow/flow/Neuman.py | 136 ++++++++++++++++++++++++++++++----------- 1 file changed, 100 insertions(+), 36 deletions(-) diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py index b3cb937..7cda613 100644 --- a/anaflow/flow/Neuman.py +++ b/anaflow/flow/Neuman.py @@ -15,52 +15,53 @@ from scipy.special import k0 from anaflow.tools.laplace import get_lap_inv from anaflow.tools.special import Shaper -from scipy.optimize import fsolve -from scipy import optimize +from scipy.optimize import root __all__ = [] -def eps_func(eps,value=0): - return eps * np.tan(eps) - value +def get_f_df(value=0): + """Get target function and its derivative.""" + if value < 0: + raise ValueError("Neuman: epsilon for root finding needs to be >0.") + def _f_df(x): + return ( + np.subtract(np.multiply(x, np.tan(x)), value), + np.tan(x) + np.divide(x, np.cos(x) ** 2), + ) -def interval_of_nth_root(n, v): - """Correct interval.""" - a, b = n * np.pi, (n + 0.5) * np.pi - 0.001 - if v < 0: # just for completeness (shift to negative branch) - a += 0.5 * np.pi - b += 0.5 * np.pi - return a, b + return _f_df def nth_root(n, v): """Get nth root of x*tan(x) = v.""" - a, b = interval_of_nth_root(n, v) - if abs(v) > 100: - return (n+0.5) * np.pi - return optimize.bisect(eps_func, a, b, (v,)) + f_df = get_f_df(v) + x0 = np.arctan(v) + n * np.pi + sol = root(f_df, x0, jac=True) + if not sol.success: + raise ValueError(f"Neuman: couldn't find {n}-th root for eps={v}") + return sol.x[0] -def neuman_unconfined_laplace( +def neuman_unconfined_partially_penetrating_laplace( s, rad, storage, transmissivity, rate, - sat_thickness=52.0, + sat_thickness=49, screen_size=11.88, well_depth=12.6, kd=0.61, specific_yield=0.26, - n_numbers=300, + n_numbers=25, ): """ - The Neuman solution. + Neuman solution for unconfined aquifers in Laplace space. - - Parameters - ---------- + Parameters + ---------- s : :class:`numpy.ndarray` Array with all "Laplace-space-points" where the function should be evaluated in the Laplace space. @@ -93,19 +94,81 @@ def neuman_unconfined_laplace( for si, se in enumerate(s): for ri, re in enumerate(rad): - if re < np.inf: - rd = re / sat_thickness - betaw = kd * (rd ** 2) - righthand = se / (((storage/specific_yield) * betaw) + se / 1e9) - roots = [nth_root(n, righthand) for n in range(n_numbers)] - for i_n in range(len(roots)): - epsilon_n = roots[i_n] - xn = (betaw * (epsilon_n ** 2) + se) ** 0.5 - res[si, ri] = (2 * k0(xn) * (np.sin(epsilon_n * (1 - (d / sat_thickness))) - np.sin(epsilon_n * - (1 - (well_depth / sat_thickness)))) * np.cos(epsilon_n * (z / sat_thickness))) / \ - (se * ((well_depth - d) / sat_thickness) * (0.5 * epsilon_n + 0.25 * np.sin(2 * epsilon_n))) - res *= rate / (2 * np.pi * transmissivity) - return res + if re == np.inf: + continue + rd = re / sat_thickness + beta = kd * (rd ** 2) + rhs = se / (((storage / specific_yield) * beta) + se / 1e9) + roots = [nth_root(n, rhs) for n in range(n_numbers)] + for eps in roots: + xn = (beta * (eps ** 2) + se) ** 0.5 + res[si, ri] += ( + 2 + * k0(xn) + * ( + np.sin(eps * (1 - (d / sat_thickness))) + - np.sin(eps * (1 - (well_depth / sat_thickness))) + ) + * np.cos(eps * (z / sat_thickness)) + ) / ( + se + * ((well_depth - d) / sat_thickness) + * (0.5 * eps + 0.25 * np.sin(2 * eps)) + ) + return res * rate / (2 * np.pi * transmissivity) + +def neuman_unconfined_fully_penetrating_laplace( + s, + rad, + storage, + transmissivity, + rate, + sat_thickness=49, + specific_yield=0.26, + n_numbers=25, +): + """ + Neuman solution for unconfined aquifers in Laplace space. + + Parameters + ---------- + s : :class:`numpy.ndarray` + Array with all "Laplace-space-points" where the function should + be evaluated in the Laplace space. + rad : :class:`numpy.ndarray` + Array with all radii where the function should be evaluated. + storage : :class:`float` + Storage of the aquifer. + transmissivity : :class:`float` + Geometric-mean transmissivity. + rate : :class:`float`, optional + Pumpingrate at the well. + sat_thickness : :class:`float`, optional + Saturated thickness of the aquifer. + kd : :class:`float`, optional + Dimensionless parameter for the conductivity. + Kz/Kr : vertical conductivity divided by horizontal conductivity + specific_yield: :class:`float`, optional + specific yield + """ + kr = transmissivity/sat_thickness + kz = kr * 0.001 + s = np.squeeze(s).reshape(-1) + rad = np.squeeze(rad).reshape(-1) + res = np.zeros(s.shape + rad.shape) + + for si, se in enumerate(s): + for ri, re in enumerate(rad): + if re == np.inf: + continue + rd = re / sat_thickness + beta = kz * (rd ** 2) / kr + rhs = se / (((storage / specific_yield) * beta) + se / ((1e9 * sat_thickness * specific_yield)/kz)) + roots = [nth_root(n, rhs) for n in range(n_numbers)] + for eps in roots: + xn = (beta * (eps ** 2) + se) ** 0.5 + res[si, ri] += (2*k0(xn) * (np.sin(eps)**2)) / (se * eps * (0.5 * eps + 0.25 * np.sin(2 * eps))) + return res * rate / (2 * np.pi * transmissivity) def neuman_unconfined( @@ -134,7 +197,7 @@ def neuman_unconfined( } kwargs.update(lap_kwargs) res = np.zeros((Input.time_no, Input.rad_no)) - lap_inv = get_lap_inv(neuman_unconfined_laplace, **kwargs) + lap_inv = get_lap_inv(neuman_unconfined_partially_penetrating_laplace, **kwargs) # call the laplace inverse function (only at time points > 0) res[Input.time_gz, :] = lap_inv(Input.time[Input.time_gz]) # reshaping results @@ -142,3 +205,4 @@ def neuman_unconfined( # add the reference head res += h_bound return res + From 728c5babe3d16a4e852140f5b64f56118680f041 Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Fri, 9 Jul 2021 15:05:43 +0200 Subject: [PATCH 7/9] Made the destinction between fully and partially penetrating wells black formatted --- anaflow/flow/Neuman.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py index 7cda613..55eaeb1 100644 --- a/anaflow/flow/Neuman.py +++ b/anaflow/flow/Neuman.py @@ -117,6 +117,7 @@ def neuman_unconfined_partially_penetrating_laplace( ) return res * rate / (2 * np.pi * transmissivity) + def neuman_unconfined_fully_penetrating_laplace( s, rad, @@ -151,7 +152,7 @@ def neuman_unconfined_fully_penetrating_laplace( specific_yield: :class:`float`, optional specific yield """ - kr = transmissivity/sat_thickness + kr = transmissivity / sat_thickness kz = kr * 0.001 s = np.squeeze(s).reshape(-1) rad = np.squeeze(rad).reshape(-1) @@ -163,11 +164,16 @@ def neuman_unconfined_fully_penetrating_laplace( continue rd = re / sat_thickness beta = kz * (rd ** 2) / kr - rhs = se / (((storage / specific_yield) * beta) + se / ((1e9 * sat_thickness * specific_yield)/kz)) + rhs = se / ( + ((storage / specific_yield) * beta) + + se / ((1e9 * sat_thickness * specific_yield) / kz) + ) roots = [nth_root(n, rhs) for n in range(n_numbers)] for eps in roots: xn = (beta * (eps ** 2) + se) ** 0.5 - res[si, ri] += (2*k0(xn) * (np.sin(eps)**2)) / (se * eps * (0.5 * eps + 0.25 * np.sin(2 * eps))) + res[si, ri] += (2 * k0(xn) * (np.sin(eps) ** 2)) / ( + se * eps * (0.5 * eps + 0.25 * np.sin(2 * eps)) + ) return res * rate / (2 * np.pi * transmissivity) @@ -197,7 +203,9 @@ def neuman_unconfined( } kwargs.update(lap_kwargs) res = np.zeros((Input.time_no, Input.rad_no)) - lap_inv = get_lap_inv(neuman_unconfined_partially_penetrating_laplace, **kwargs) + lap_inv = get_lap_inv( + neuman_unconfined_partially_penetrating_laplace, **kwargs + ) # call the laplace inverse function (only at time points > 0) res[Input.time_gz, :] = lap_inv(Input.time[Input.time_gz]) # reshaping results @@ -205,4 +213,3 @@ def neuman_unconfined( # add the reference head res += h_bound return res - From 1a3f275f992c8e4a24228284578e7298bcf0d3a9 Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Tue, 10 Aug 2021 10:31:05 +0200 Subject: [PATCH 8/9] Added the latest rootfinding and corrected both the inint files --- anaflow/__init__.py | 4 ++++ anaflow/flow/Neuman.py | 36 +++++++++++++++++++----------------- anaflow/flow/__init__.py | 8 +++++++- 3 files changed, 30 insertions(+), 18 deletions(-) diff --git a/anaflow/__init__.py b/anaflow/__init__.py index 67ad3c5..62b7946 100644 --- a/anaflow/__init__.py +++ b/anaflow/__init__.py @@ -99,6 +99,8 @@ neuman2004_steady, ext_grf, ext_grf_steady, + neuman_unconfined_partially_penetrating_laplace, + neuman_unconfined_fully_penetrating_laplace, ) from anaflow.tools import ( get_lap_inv, @@ -136,4 +138,6 @@ "step_f", "specialrange", "specialrange_cut", + "neuman_unconfined_partially_penetrating_laplace", + "neuman_unconfined_fully_penetrating_laplace", ] diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py index 55eaeb1..412e409 100644 --- a/anaflow/flow/Neuman.py +++ b/anaflow/flow/Neuman.py @@ -7,8 +7,11 @@ The following functions are provided .. autosummary:: - neuman_unconfined_laplace - neuman_unconfined + get_f_df + nth_root + neuman_unconfined_partially_penetrating_laplace + neuman_unconfined_fully_penetrating_laplace + neuman_unconfined """ # pylint: disable=C0103 import numpy as np @@ -16,33 +19,32 @@ from anaflow.tools.laplace import get_lap_inv from anaflow.tools.special import Shaper from scipy.optimize import root +from scipy.optimize import root_scalar __all__ = [] - -def get_f_df(value=0): - """Get target function and its derivative.""" +def get_f_df_df2(value=0): + """Get target function and its first two derivatives.""" if value < 0: raise ValueError("Neuman: epsilon for root finding needs to be >0.") - def _f_df(x): + def _f_df_df2(x): return ( np.subtract(np.multiply(x, np.tan(x)), value), np.tan(x) + np.divide(x, np.cos(x) ** 2), + 2 * (np.multiply(x, np.tan(x)) + 1.0) * np.cos(x) ** -2, ) - return _f_df - + return _f_df_df2 def nth_root(n, v): - """Get nth root of x*tan(x) = v.""" - f_df = get_f_df(v) - x0 = np.arctan(v) + n * np.pi - sol = root(f_df, x0, jac=True) - if not sol.success: + """Get n-th root of x*tan(x) = v.""" + x0 = np.sqrt(v) if (v < 1 and n < 1) else np.arctan(v) + n * np.pi + f = get_f_df_df2(v) + sol = root_scalar(f, x0=x0, fprime=True, fprime2=True) + if not sol.converged: raise ValueError(f"Neuman: couldn't find {n}-th root for eps={v}") - return sol.x[0] - + return sol.root def neuman_unconfined_partially_penetrating_laplace( s, @@ -84,7 +86,7 @@ def neuman_unconfined_partially_penetrating_laplace( Dimensionless parameter for the conductivity. Kz/Kr : vertical conductivity divided by horizontal conductivity specific_yield: :class:`float`, optional - specific yield + Specific yield """ z = sat_thickness - well_depth d = well_depth - screen_size @@ -150,7 +152,7 @@ def neuman_unconfined_fully_penetrating_laplace( Dimensionless parameter for the conductivity. Kz/Kr : vertical conductivity divided by horizontal conductivity specific_yield: :class:`float`, optional - specific yield + Specific yield """ kr = transmissivity / sat_thickness kz = kr * 0.001 diff --git a/anaflow/flow/__init__.py b/anaflow/flow/__init__.py index 7364780..83f597c 100644 --- a/anaflow/flow/__init__.py +++ b/anaflow/flow/__init__.py @@ -72,11 +72,15 @@ neuman2004_steady, ) from anaflow.flow.ext_grf_model import ext_grf, ext_grf_steady +from anaflow.flow.Neuman import ( + neuman_unconfined_partially_penetrating_laplace, + neuman_unconfined_fully_penetrating_laplace, +) __all__ = [ "thiem", "theis", - "grf_model", + "grf", "ext_thiem_2d", "ext_thiem_3d", "ext_thiem_tpl", @@ -90,4 +94,6 @@ "grf", "ext_grf", "ext_grf_steady", + "neuman_unconfined_partially_penetrating_laplace", + "neuman_unconfined_fully_penetrating_laplace", ] From 6db00fccabbcd2fc35710cdcb981c60b756ae84b Mon Sep 17 00:00:00 2001 From: JarnoHerr Date: Tue, 10 Aug 2021 10:36:51 +0200 Subject: [PATCH 9/9] Added the latest rootfinding and corrected both the init files --- anaflow/__init__.py | 6 ++---- anaflow/flow/Neuman.py | 1 - anaflow/flow/__init__.py | 9 +++------ 3 files changed, 5 insertions(+), 11 deletions(-) diff --git a/anaflow/__init__.py b/anaflow/__init__.py index 62b7946..d98f84f 100644 --- a/anaflow/__init__.py +++ b/anaflow/__init__.py @@ -99,8 +99,7 @@ neuman2004_steady, ext_grf, ext_grf_steady, - neuman_unconfined_partially_penetrating_laplace, - neuman_unconfined_fully_penetrating_laplace, + neuman_unconfined, ) from anaflow.tools import ( get_lap_inv, @@ -138,6 +137,5 @@ "step_f", "specialrange", "specialrange_cut", - "neuman_unconfined_partially_penetrating_laplace", - "neuman_unconfined_fully_penetrating_laplace", + "neuman_unconfined", ] diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py index 412e409..a5499c1 100644 --- a/anaflow/flow/Neuman.py +++ b/anaflow/flow/Neuman.py @@ -18,7 +18,6 @@ from scipy.special import k0 from anaflow.tools.laplace import get_lap_inv from anaflow.tools.special import Shaper -from scipy.optimize import root from scipy.optimize import root_scalar __all__ = [] diff --git a/anaflow/flow/__init__.py b/anaflow/flow/__init__.py index 83f597c..65e3a23 100644 --- a/anaflow/flow/__init__.py +++ b/anaflow/flow/__init__.py @@ -72,10 +72,8 @@ neuman2004_steady, ) from anaflow.flow.ext_grf_model import ext_grf, ext_grf_steady -from anaflow.flow.Neuman import ( - neuman_unconfined_partially_penetrating_laplace, - neuman_unconfined_fully_penetrating_laplace, -) +from anaflow.flow.Neuman import neuman_unconfined + __all__ = [ "thiem", @@ -94,6 +92,5 @@ "grf", "ext_grf", "ext_grf_steady", - "neuman_unconfined_partially_penetrating_laplace", - "neuman_unconfined_fully_penetrating_laplace", + "neuman_unconfined" ]