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/* 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 diff --git a/anaflow/__init__.py b/anaflow/__init__.py index 67ad3c5..d98f84f 100644 --- a/anaflow/__init__.py +++ b/anaflow/__init__.py @@ -99,6 +99,7 @@ neuman2004_steady, ext_grf, ext_grf_steady, + neuman_unconfined, ) from anaflow.tools import ( get_lap_inv, @@ -136,4 +137,5 @@ "step_f", "specialrange", "specialrange_cut", + "neuman_unconfined", ] diff --git a/anaflow/flow/Neuman.py b/anaflow/flow/Neuman.py new file mode 100644 index 0000000..a5499c1 --- /dev/null +++ b/anaflow/flow/Neuman.py @@ -0,0 +1,216 @@ +# -*- coding: utf-8 -*- +""" +Anaflow subpackage providing the Neuman equation for homogeneous aquifer. + +.. currentmodule:: anaflow.flow.Neuman + +The following functions are provided + +.. autosummary:: + get_f_df + nth_root + neuman_unconfined_partially_penetrating_laplace + neuman_unconfined_fully_penetrating_laplace + neuman_unconfined +""" +# pylint: disable=C0103 +import numpy as np +from scipy.special import k0 +from anaflow.tools.laplace import get_lap_inv +from anaflow.tools.special import Shaper +from scipy.optimize import root_scalar + +__all__ = [] + +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_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_df2 + +def nth_root(n, v): + """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.root + +def neuman_unconfined_partially_penetrating_laplace( + s, + rad, + storage, + transmissivity, + rate, + sat_thickness=49, + screen_size=11.88, + well_depth=12.6, + kd=0.61, + 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. + sat_thickness : :class:`float`, optional + Saturated thickness of the aquifer. + rate : :class:`float`, optional + Pumpingrate at the well. Default: -1e-4 + 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 + 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 = 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) + + for si, se in enumerate(s): + for ri, re in enumerate(rad): + 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( + 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) + 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_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 + res = Input.reshape(res) + # add the reference head + res += h_bound + return res diff --git a/anaflow/flow/__init__.py b/anaflow/flow/__init__.py index 7364780..65e3a23 100644 --- a/anaflow/flow/__init__.py +++ b/anaflow/flow/__init__.py @@ -72,11 +72,13 @@ neuman2004_steady, ) from anaflow.flow.ext_grf_model import ext_grf, ext_grf_steady +from anaflow.flow.Neuman import neuman_unconfined + __all__ = [ "thiem", "theis", - "grf_model", + "grf", "ext_thiem_2d", "ext_thiem_3d", "ext_thiem_tpl", @@ -90,4 +92,5 @@ "grf", "ext_grf", "ext_grf_steady", + "neuman_unconfined" ]