diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 3f18d32..ad90c47 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,34 +12,33 @@ on: # Allows you to run this workflow manually from the Actions tab workflow_dispatch: -env: - # needed by coveralls - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - jobs: source_check: name: source check runs-on: ubuntu-latest - strategy: - fail-fast: false steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - - name: Set up Python 3.9 - uses: actions/setup-python@v2 + - name: Set up Python 3.11 + uses: actions/setup-python@v5 with: - python-version: 3.9 + python-version: 3.11 - name: Install dependencies run: | python -m pip install --upgrade pip - pip install --editable .[check] + pip install --editable .[check,test] + pip install "coveralls>=3.0.0" - name: black check run: | python -m black --check --diff --color . + - name: black preview + run: | + python -m black --preview --diff --color . + - name: isort check run: | python -m isort --check --diff --color . @@ -48,53 +47,60 @@ jobs: run: | python -m pylint src/anaflow/ + - name: coveralls check + run: | + python -m pytest --cov anaflow --cov-report term-missing -v tests/ + python -m coveralls --service=github + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + 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.7', '3.8', '3.9', '3.10', '3.11'] + os: [ubuntu-latest, windows-latest, macos-13] + python-version: ['3.8', '3.9', '3.10', '3.11', '3.12', '3.13'] steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 with: fetch-depth: '0' - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v2 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} - name: Install dependencies run: | python -m pip install --upgrade pip - pip install build coveralls>=3.0.0 + pip install build pip install --editable .[test] - name: Run tests run: | python -m pytest --cov anaflow --cov-report term-missing -v tests/ - python -m coveralls --service=github - name: Build sdist run: | python -m build - - uses: actions/upload-artifact@v2 - if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9' + - uses: actions/upload-artifact@v4 + if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.11' with: - path: dist + name: dist + path: dist/* upload_to_pypi: needs: [build_sdist] runs-on: ubuntu-latest steps: - - uses: actions/download-artifact@v2 + - uses: actions/download-artifact@v4 with: - name: artifact + name: dist path: dist - name: Publish to Test PyPI diff --git a/.gitignore b/.gitignore index ce93941..9b5c028 100644 --- a/.gitignore +++ b/.gitignore @@ -63,6 +63,7 @@ instance/ # Sphinx documentation docs/_build/ +docs/source/sg_execution_times.rst # PyBuilder target/ diff --git a/.readthedocs.yml b/.readthedocs.yml index f88c11a..4291bf8 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -1,12 +1,16 @@ version: 2 +build: + os: ubuntu-22.04 + tools: + python: "3.11" + sphinx: configuration: docs/source/conf.py -formats: all +formats: [pdf] python: - version: 3.7 install: - method: pip path: . diff --git a/CHANGELOG.md b/CHANGELOG.md index 4562503..4880da7 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,24 @@ All notable changes to **AnaFlow** will be documented in this file. +## [1.2.0] - 2025-05 + +See [#12](https://github.com/GeoStat-Framework/AnaFlow/pull/12) + +### Enhancements +- added solutions based on the effective transmissivity for the ["Integral" variogram model](https://geostat-framework.readthedocs.io/projects/gstools/en/v1.7.0/api/gstools.covmodel.Integral.html): + - `ext_thiem_int`: steady state solution + - `ext_thiem_int_3d`: steady state solution incorporating vertical anisotropy + - `ext_theis_int`: transient solution + - `ext_theis_int_3d`: transient solution incorporating vertical anisotropy +- added `fix_T_well` and `fix_K_well` bool flag to transient heterogeneous solutions to be able to set if the well value for the effective transmissivity/conductivity should be determined from the limit (`True`) or from the upscaled value in the first ring segment (`False`, default) + - **breaking**: the previous behavior was effectively this set to `True`, which for steep effective curves resulted in an increasing error in the effective head near the well + +### Changes +- updated docs (use myst parser for markdown files, only generate html and pdf) +- updated CI (fixed artifacts up-/download action, see #14) +- use hatchling as build backend + ## [1.1.0] - 2023-04 @@ -88,6 +106,7 @@ Containing: - lap_transgwflow_cyl - Solution for a diskmodel in laplace inversion +[1.2.0]: https://github.com/GeoStat-Framework/AnaFlow/compare/v1.1.0...v1.2.0 [1.1.0]: https://github.com/GeoStat-Framework/AnaFlow/compare/v1.0.1...v1.1.0 [1.0.1]: https://github.com/GeoStat-Framework/AnaFlow/compare/v1.0.0...v1.0.1 [1.0.0]: https://github.com/GeoStat-Framework/AnaFlow/compare/v0.4.0...v1.0.0 diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100755 index 15fa4ed..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,4 +0,0 @@ -prune ** -recursive-include tests *.py -recursive-include src/anaflow *.py -include LICENSE README.md pyproject.toml diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index ab37940..58ac813 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -1 +1,2 @@ -.. mdinclude:: ../../CHANGELOG.md +.. include:: ../../CHANGELOG.md + :parser: myst_parser.docutils_ \ No newline at end of file diff --git a/docs/source/conf.py b/docs/source/conf.py index f638bc5..a455e61 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -55,7 +55,7 @@ def setup(app): "sphinx.ext.napoleon", # parameters look better than with numpydoc only "numpydoc", "sphinx_gallery.gen_gallery", - "m2r2", + "myst_parser", ] # autosummaries from source-files @@ -87,8 +87,12 @@ def setup(app): # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: +source_suffix = { + ".rst": "restructuredtext", + ".md": "markdown", +} # source_suffix = ['.rst', '.md'] -source_suffix = ".rst" +# source_suffix = ".rst" # The master toctree document. # --> this is the sitemap (or content-list in latex -> needs a heading) @@ -141,7 +145,7 @@ def setup(app): # 'canonical_url': '', # 'analytics_id': '', "logo_only": False, - "display_version": True, + "version_selector": True, "prev_next_buttons_location": "top", # 'style_external_links': False, # 'vcs_pageview_mode': '', @@ -267,7 +271,7 @@ def setup(app): # Remove the "Download all examples" button from the top level gallery "download_all_examples": False, # Sort gallery example by file name instead of number of lines (default) - "within_subsection_order": FileNameSortKey, + "within_subsection_order": "FileNameSortKey", # directory where function granular galleries are stored "backreferences_dir": None, # Modules for which function level galleries are created. In diff --git a/examples/01_call_theis.py b/examples/01_call_theis.py index d7281dc..6bc4a3c 100644 --- a/examples/01_call_theis.py +++ b/examples/01_call_theis.py @@ -7,6 +7,7 @@ Reference: `Theis 1935 `__ """ + import numpy as np from matplotlib import pyplot as plt diff --git a/examples/02_call_ext_theis2d.py b/examples/02_call_ext_theis2d.py index c37490d..769de8b 100644 --- a/examples/02_call_ext_theis2d.py +++ b/examples/02_call_ext_theis2d.py @@ -12,6 +12,7 @@ Reference: `Zech et. al. 2016 `__ """ + import numpy as np from matplotlib import pyplot as plt @@ -51,12 +52,8 @@ label_TG = "Theis($T_G$)" if i == 0 else None label_TH = "Theis($T_H$)" if i == 0 else None label_ef = "extended Theis" if i == 0 else None - plt.plot( - rad, head_TG[i], label=label_TG, color="C" + str(i), linestyle="--" - ) - plt.plot( - rad, head_TH[i], label=label_TH, color="C" + str(i), linestyle=":" - ) + plt.plot(rad, head_TG[i], label=label_TG, color="C" + str(i), linestyle="--") + plt.plot(rad, head_TH[i], label=label_TH, color="C" + str(i), linestyle=":") plt.plot(rad, head_ef[i], label=label_ef, color="C" + str(i)) time_ticks.append(head_ef[i][-1]) diff --git a/examples/03_call_ext_theis3d.py b/examples/03_call_ext_theis3d.py index 19d56bc..af91897 100644 --- a/examples/03_call_ext_theis3d.py +++ b/examples/03_call_ext_theis3d.py @@ -14,6 +14,7 @@ Reference: `Müller 2015 `__ """ + import numpy as np from matplotlib import pyplot as plt @@ -59,12 +60,8 @@ label_TG = "Theis($K_{efu}$)" if i == 0 else None label_TH = "Theis($K_H$)" if i == 0 else None label_ef = "extended Theis 3D" if i == 0 else None - plt.plot( - rad, head_Kefu[i], label=label_TG, color="C" + str(i), linestyle="--" - ) - plt.plot( - rad, head_KH[i], label=label_TH, color="C" + str(i), linestyle=":" - ) + plt.plot(rad, head_Kefu[i], label=label_TG, color="C" + str(i), linestyle="--") + plt.plot(rad, head_KH[i], label=label_TH, color="C" + str(i), linestyle=":") plt.plot(rad, head_ef[i], label=label_ef, color="C" + str(i)) time_ticks.append(head_ef[i][-1]) diff --git a/examples/04_call_ext_theis_tpl.py b/examples/04_call_ext_theis_tpl.py index 82271b0..e234d8a 100755 --- a/examples/04_call_ext_theis_tpl.py +++ b/examples/04_call_ext_theis_tpl.py @@ -13,6 +13,7 @@ Reference: (not yet published) """ + import numpy as np from matplotlib import pyplot as plt @@ -62,12 +63,8 @@ label_TG = "Theis($K_G$)" if i == 0 else None label_TH = "Theis($K_H$)" if i == 0 else None label_ef = "extended Theis TPL 2D" if i == 0 else None - plt.plot( - rad, head_KG[i], label=label_TG, color="C" + str(i), linestyle="--" - ) - plt.plot( - rad, head_KH[i], label=label_TH, color="C" + str(i), linestyle=":" - ) + plt.plot(rad, head_KG[i], label=label_TG, color="C" + str(i), linestyle="--") + plt.plot(rad, head_KH[i], label=label_TH, color="C" + str(i), linestyle=":") plt.plot(rad, head_ef[i], label=label_ef, color="C" + str(i)) time_ticks.append(head_ef[i][-1]) diff --git a/examples/05_call_neuman2004.py b/examples/05_call_neuman2004.py index c35c6cb..02c7de7 100755 --- a/examples/05_call_neuman2004.py +++ b/examples/05_call_neuman2004.py @@ -14,6 +14,7 @@ Reference: `Neuman 2004 `__ """ + import numpy as np from matplotlib import pyplot as plt @@ -61,12 +62,8 @@ label_TG = "Theis($T_G$)" if i == 0 else None label_TH = "Theis($T_H$)" if i == 0 else None label_ef = "transient Neuman [2004]" if i == 0 else None - plt.plot( - rad, head_TG[i], label=label_TG, color="C" + str(i), linestyle="--" - ) - plt.plot( - rad, head_TH[i], label=label_TH, color="C" + str(i), linestyle=":" - ) + plt.plot(rad, head_TG[i], label=label_TG, color="C" + str(i), linestyle="--") + plt.plot(rad, head_TH[i], label=label_TH, color="C" + str(i), linestyle=":") plt.plot(rad, head_ef[i], label=label_ef, color="C" + str(i)) time_ticks.append(head_ef[i][-1]) diff --git a/examples/06_compare_extthiem2d_grfsteady.py b/examples/06_compare_extthiem2d_grfsteady.py index 2540e81..1edae14 100755 --- a/examples/06_compare_extthiem2d_grfsteady.py +++ b/examples/06_compare_extthiem2d_grfsteady.py @@ -11,6 +11,7 @@ - `Schneider & Attinger 2008 `__ - `Zech & Attinger 2016 `__ """ + import numpy as np from matplotlib import pyplot as plt diff --git a/examples/07_compare_extthiem3d_grfsteady.py b/examples/07_compare_extthiem3d_grfsteady.py index 208c273..aa3b4e0 100755 --- a/examples/07_compare_extthiem3d_grfsteady.py +++ b/examples/07_compare_extthiem3d_grfsteady.py @@ -8,6 +8,7 @@ Reference: `Zech et. al. 2012 `__ """ + import numpy as np from matplotlib import pyplot as plt diff --git a/examples/08_compare_extthiem2d_neuman.py b/examples/08_compare_extthiem2d_neuman.py index 4e442c8..67c8da7 100755 --- a/examples/08_compare_extthiem2d_neuman.py +++ b/examples/08_compare_extthiem2d_neuman.py @@ -13,6 +13,7 @@ - `Neuman 2004 `__ - `Zech & Attinger 2016 `__ """ + import numpy as np from matplotlib import pyplot as plt diff --git a/examples/09_compare_exttheis2d_neuman.py b/examples/09_compare_exttheis2d_neuman.py index 25fbd4e..c388a82 100755 --- a/examples/09_compare_exttheis2d_neuman.py +++ b/examples/09_compare_exttheis2d_neuman.py @@ -13,6 +13,7 @@ - `Neuman 2004 `__ - `Zech et. al. 2016 `__ """ + import numpy as np from matplotlib import pyplot as plt @@ -40,9 +41,7 @@ plt.plot(rad, head2[i], label=label2, color="C" + str(i), linestyle="--") time_ticks.append(head1[i][-1]) -plt.title( - "$T_G={}$, $\sigma^2={}$, $\ell={}$, $S={}$".format(TG, var, len_scale, S) -) +plt.title(r"$T_G={}$, $\sigma^2={}$, $\ell={}$, $S={}$".format(TG, var, len_scale, S)) plt.xlabel("r in [m]") plt.ylabel("h in [m]") plt.legend() diff --git a/examples/10_convergence_ext_theis_tpl.py b/examples/10_convergence_ext_theis_tpl.py index dee18fd..25138fe 100755 --- a/examples/10_convergence_ext_theis_tpl.py +++ b/examples/10_convergence_ext_theis_tpl.py @@ -7,6 +7,7 @@ Reference: (not yet published) """ + import numpy as np from matplotlib import pyplot as plt @@ -25,9 +26,7 @@ S = 1e-4 # storativity rate = -1e-4 # pumping rate -head1 = ext_thiem_tpl( - rad, r_ref, KG, len_scale, hurst, var, dim=dim, rate=rate -) +head1 = ext_thiem_tpl(rad, r_ref, KG, len_scale, hurst, var, dim=dim, rate=rate) head2 = ext_theis_tpl( time, rad, S, KG, len_scale, hurst, var, dim=dim, rate=rate, r_bound=r_ref ) diff --git a/examples/11_convergence_ext_grf.py b/examples/11_convergence_ext_grf.py index 518aed7..ad20f4e 100755 --- a/examples/11_convergence_ext_grf.py +++ b/examples/11_convergence_ext_grf.py @@ -9,6 +9,7 @@ Reference: `Barker 1988 `__ """ + import numpy as np from matplotlib import pyplot as plt @@ -30,9 +31,7 @@ plt.plot(rad, head1, label="Ext GRF steady") plt.plot(rad, head2, label="Ext GRF (t={})".format(time), linestyle="--") -plt.plot( - rad, head3, label="GRF quasi-steady (t={})".format(time), linestyle=":" -) +plt.plot(rad, head3, label="GRF quasi-steady (t={})".format(time), linestyle=":") plt.xlabel("r in [m]") plt.ylabel("h in [m]") diff --git a/examples/12_compare_theis_quasi_steady.py b/examples/12_compare_theis_quasi_steady.py index ea89012..81cf687 100755 --- a/examples/12_compare_theis_quasi_steady.py +++ b/examples/12_compare_theis_quasi_steady.py @@ -5,6 +5,7 @@ The quasi steady is reached, when the radial shape of the drawdown in not changing anymore. """ + import numpy as np from matplotlib import pyplot as plt @@ -21,12 +22,8 @@ transmissivity=1e-4, rate=-1e-4, ) -head1 = ( - theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4) - head_ref -) -head2 = theis( - time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4, r_bound=r_ref -) +head1 = theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4) - head_ref +head2 = theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4, r_bound=r_ref) head3 = thiem(rad, r_ref, transmissivity=1e-4, rate=-1e-4) for i, step in enumerate(time): diff --git a/examples/13_self_defined_transmissivity.py b/examples/13_self_defined_transmissivity.py index 72fe491..71899c0 100755 --- a/examples/13_self_defined_transmissivity.py +++ b/examples/13_self_defined_transmissivity.py @@ -15,6 +15,7 @@ Reference: (not yet published) """ + import matplotlib.gridspec as gridspec import numpy as np from matplotlib import pyplot as plt @@ -78,9 +79,7 @@ def cond(rad, K_far, K_well, len_scale): rad_lin = np.linspace(rad[0], rad[-1], 1000) ax1.plot(rad_lin, step_f(rad_lin, R_part, K_part), label="step Conductivity") -ax1.plot( - rad_lin, cond(rad_lin, K_far, K_well, len_scale), label="Conductivity" -) +ax1.plot(rad_lin, cond(rad_lin, K_far, K_well, len_scale), label="Conductivity") ax1.set_yticks([K_well, K_far]) ax1.set_ylabel(r"$K$ in $[\frac{m}{s}]$") plt.setp(ax1.get_xticklabels(), visible=False) diff --git a/examples/14_interval_theis.py b/examples/14_interval_theis.py index ad9be56..084765d 100755 --- a/examples/14_interval_theis.py +++ b/examples/14_interval_theis.py @@ -8,6 +8,7 @@ Unfortunatly the Stehfest algorithm is not suitable for this kind of solution, which is demonstrated in the following script. """ + import numpy as np from matplotlib import pyplot as plt diff --git a/examples/15_accruing_theis.py b/examples/15_accruing_theis.py index 523d041..d285d20 100755 --- a/examples/15_accruing_theis.py +++ b/examples/15_accruing_theis.py @@ -7,6 +7,7 @@ This could be interpreted as that the water pump needs a certain time to reach its constant rate state. """ + import matplotlib.gridspec as gridspec import numpy as np from matplotlib import pyplot as plt diff --git a/examples/16_integral_model.py b/examples/16_integral_model.py new file mode 100644 index 0000000..a28dd31 --- /dev/null +++ b/examples/16_integral_model.py @@ -0,0 +1,56 @@ +r""" +The integral variogram model +============================ +""" + +import numpy as np +from matplotlib import pyplot as plt + +import anaflow as af + + +def dashes(i=1, max_n=12, width=1): + return i * [width, width] + [max_n * 2 * width - 2 * i * width, width] + + +rad = np.linspace(0, 3, 1000) +rough = [0.01, 0.1, 0.25, 0.5, 1.0, 1.5, 2.0, 10.0] +# rescale for integral scale of 1 +length = [1 / (nu * np.sqrt(np.pi) / (2 * nu + 2.0)) for nu in rough] +ln_gau = 2 / np.sqrt(np.pi) +var = 1 +TG = 1e-4 +TH = TG * np.exp(-var / 2) + +for i, (ln, nu) in enumerate(zip(length, rough)): + t = 1e4 * af.tools.Int_CG(rad, TG, var, ln, nu) + plt.plot( + rad, + t, + label=r"$T^{Int}_{CG}$($\nu$" + f"={nu:.4})", + dashes=dashes(i), + color="C0" if nu < 1 else "C2", + ) +plt.plot( + rad, + af.tools.T_CG(rad, TG, var, ln_gau) / TG, + label=r"$T_{CG}$ - Gaussian", + color="C1", +) +plt.title(f"$T_G=1$ cm²/s, $\\sigma^2={var}$") +plt.xlabel(r"r / $\ell$") +plt.ylabel(r"T / cm²s⁻¹") +plt.grid() +ax = plt.gca() +ax.legend() +ax.ticklabel_format(axis="y", style="scientific", scilimits=[-3, 0], useMathText=True) +ax.locator_params(tight=True, nbins=5) + +ylim = ax.get_ylim() +ax2 = ax.twinx() +ax2.set_ylim(ylim) +ax2.set_yticks([TH * 1e4, TG * 1e4]) +ax2.set_yticklabels([r"$T_H$", r"$T_G$"]) + +plt.tight_layout() +plt.show() diff --git a/examples/17_compare_nugget.py b/examples/17_compare_nugget.py new file mode 100644 index 0000000..764d0e5 --- /dev/null +++ b/examples/17_compare_nugget.py @@ -0,0 +1,86 @@ +""" +Impact of roughness on effective drawdown on anti-persistent fields +=================================================================== + +When the field roughness approaches zero, the drawdown approaches +the homogeneous solution for the geometric mean transmissivity. +""" + +import numpy as np +from matplotlib import pyplot as plt + +import anaflow as af + +min_s = 60 +hour_s = min_s * 60 +day_s = hour_s * 24 + + +def dashes(i=1, max_n=12, width=1): + return i * [width, width] + [max_n * 2 * width - 2 * i * width, width] + + +time_labels = ["10 s", "30 min", "steady\nshape"] +time = [10, 30 * min_s, 10 * day_s] # 10s, 30min, 10d +rad = np.geomspace(0.05, 4) # radius from the pumping well in [0, 4] + +TG = 1e-4 # the geometric mean of the transmissivity +var = 2.25 # variance of the log-transmissivity +len_scale = 10 # correlation length of the log-transmissivity + +# different roughness values from very rough to very smooth +rough = [0.01, 0.1, 0.25, 0.5, 1.0] +S = 1e-4 # storativity +rate = -1e-4 # pumping rate +heads = [] +hom_theis = af.theis(time, rad, S, TG, rate=rate) +ext_theis = af.ext_theis_2d( + time, rad, S, TG, var, len_scale / np.sqrt(np.pi) * 2, rate=rate +) + +for roughness in rough: + # rescale to constant integral scale + ln = len_scale / (roughness * np.sqrt(np.pi) / (2 * roughness + 2.0)) + # ln = len_scale + heads.append( + af.ext_theis_int(time, rad, S, TG, var, ln, roughness, rate=rate, parts=50) + ) + +time_ticks = [] +for i, step in enumerate(time): + for j, roughness in enumerate(rough): + label = ( + r"Theis$^{Int}_{CG}(\nu=" + f"{roughness:.4})$" + if i == len(time) - 1 + else None + ) + plt.plot( + rad, + heads[j][i], + label=label, + color=f"C0", + dashes=dashes(j), + alpha=(1 + i) / len(time), + ) + time_ticks.append(hom_theis[i][-1]) + label = "Theis($T_G$)" if i == len(time) - 1 else None + plt.plot(rad, hom_theis[i], label=label, color="k") # , dashes=dashes(1)) + + +plt.title( + f"$T_G=1$ cm²/s, $\\sigma^2={var}$, $\\ell={len_scale}$ m, $S={S}$, $Q={-rate*1000}$ L/s" +) +plt.xlabel("r / m") +plt.ylabel("h / m") +plt.grid() +plt.legend() +plt.gca().locator_params(tight=True, nbins=5) +plt.gca().set_ylim([-3.2, 0.1]) +plt.gca().set_xlim([0, rad[-1]]) +ylim = plt.gca().get_ylim() +ax2 = plt.gca().twinx() +ax2.set_yticks(time_ticks) +ax2.set_yticklabels(time_labels) +ax2.set_ylim(ylim) +plt.tight_layout() +plt.show() diff --git a/examples/18_compare_roughness.py b/examples/18_compare_roughness.py new file mode 100644 index 0000000..4b66d8b --- /dev/null +++ b/examples/18_compare_roughness.py @@ -0,0 +1,87 @@ +r""" +Impact of roughness on effective drawdown on persistent fields +============================================================== + +When the field gets smoother, the drawdown approaches +the effective theis solution for a gaussian correlation structure. +""" + +import numpy as np +from matplotlib import pyplot as plt + +import anaflow as af + +min_s = 60 +hour_s = min_s * 60 +day_s = hour_s * 24 + + +def dashes(i=1, max_n=12, width=1): + return i * 2 * [width] + [max_n * 2 * width - 2 * i * width, width] + + +time_labels = ["10 s", "30 min", "steady\nshape"] +time = [10, 30 * min_s, 10 * day_s] # 10s, 30min, 10days + +rad = np.geomspace(0.05, 4) # radius from the pumping well in [0, 4] + +TG = 1e-4 # the geometric mean of the transmissivity +var = 2.25 # 3 # correlation length of the log-transmissivity +len_scale = 10 # variance of the log-transmissivity + +# different roughness values from very rough to very smooth +rough = [1.0, 1.5, 2.0, 10.0] + +S = 1e-4 # storativity +rate = -1e-4 # pumping rate +heads = [] +hom_theis = af.theis(time, rad, S, TG, rate=rate) +ext_theis = af.ext_theis_2d( + time, rad, S, TG, var, len_scale / np.sqrt(np.pi) * 2, rate=rate +) + +for roughness in rough: + # rescale to constant integral scale + ln = len_scale / (roughness * np.sqrt(np.pi) / (2 * roughness + 2.0)) + heads.append( + af.ext_theis_int(time, rad, S, TG, var, ln, roughness, rate=rate, parts=50) + ) + +time_ticks = [] +for i, step in enumerate(time): + for j, roughness in enumerate(rough): + label = ( + r"Theis$^{Int}_{CG}(\nu=" + f"{roughness:.4})$" + if i == len(time) - 1 + else None + ) + plt.plot( + rad, + heads[j][i], + label=label, + color=f"C2", + dashes=dashes(j), + alpha=(1 + i) / len(time), + ) + time_ticks.append(heads[-1][i][-1]) + label = "extended Theis" if i == len(time) - 1 else None + plt.plot(rad, ext_theis[i], label=label, color="k") + + +plt.title( + f"$T_G=1$ cm²/s, $\\sigma^2={var}$, $\\ell={len_scale}$ m, $S={S}$, $Q={-rate*1000}$ L/s" +) +plt.xlabel("r / m") +plt.ylabel("h / m") +plt.grid() +plt.legend() +plt.gca().locator_params(tight=True, nbins=5) +plt.gca().set_ylim([-3.2, 0.1]) +plt.gca().set_xlim([0, rad[-1]]) +ylim = plt.gca().get_ylim() +ax2 = plt.gca().twinx() +ax2.set_yticks(time_ticks) +ax2.set_yticklabels(time_labels) +ax2.set_ylim(ylim) +plt.tight_layout() +plt.show() diff --git a/examples/19_convergence_ext_theis_int.py b/examples/19_convergence_ext_theis_int.py new file mode 100644 index 0000000..07d2891 --- /dev/null +++ b/examples/19_convergence_ext_theis_int.py @@ -0,0 +1,62 @@ +r""" +Convergence of the extended Theis solutions for the Integral model +================================================================== + +Here we set an outer boundary to the transient solution, so this condition +coincides with the references head of the steady solution. +""" + +import numpy as np +from matplotlib import pyplot as plt + +from anaflow import ext_theis_int, ext_thiem_int + +time = 86400 # time point for steady state (1 day) +rad = np.geomspace(0.05, 5) / 5 # radius from the pumping well in [0, 1] +r_ref = 2 # reference radius + +KG = 1e-4 # the geometric mean of the transmissivity +len_scale = 1 # correlation length of the log-transmissivity +roughness = 1 # roughness parameter +var = 1 # variance of the log-transmissivity +S = 1e-4 # storativity +rate = -1e-4 # pumping rate +dim = 2 + +head1 = ext_thiem_int(rad, r_ref, KG, var, len_scale, roughness, dim=dim, rate=rate) +head2 = ext_theis_int( + time, + rad, + S, + KG, + var, + len_scale, + roughness, + dim=dim, + rate=rate, + r_bound=r_ref, +) + +plt.plot(rad, head1, label=r"Thiem$^{Int}_{CG}$", linewidth=3, color="k") +plt.plot( + rad, + head2, + label=r"Theis$^{Int}_{CG}$" + f"(t=1 day)", + linestyle="--", + linewidth=3, + color="C1", +) +plt.title( + f"$T_G=1$ cm²/s, $\\sigma^2={var}$, $\\nu={roughness}$, $S={S}$, $Q=0.1$ L/s, " + + r"$r_{ref}=2\ell$" +) + +plt.xlabel(r"r / $\ell$") +plt.ylabel("h / m") +plt.grid() +plt.gca().set_ylim([-1.25, 0]) +plt.gca().set_xlim([0, 1]) +plt.gca().locator_params(tight=True, nbins=5) +plt.tight_layout() +plt.legend(loc="lower right") +plt.show() diff --git a/pyproject.toml b/pyproject.toml index cb46f56..1f638a0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,16 +1,16 @@ [build-system] requires = [ - "setuptools>=64", - "setuptools_scm>=7", + "hatchling>=1.8.0", + "hatch-vcs", ] -build-backend = "setuptools.build_meta" +build-backend = "hatchling.build" [project] -requires-python = ">=3.7" +requires-python = ">=3.8" name = "anaflow" authors = [{name = "Sebastian Mueller", email = "sebastian.mueller@ufz.de"}] readme = "README.md" -license = {text = "MIT"} +license = "MIT" dynamic = ["version"] description = "AnaFlow - analytical solutions for the groundwater-flow equation." classifiers = [ @@ -29,34 +29,36 @@ classifiers = [ "Programming Language :: Python", "Programming Language :: Python :: 3", "Programming Language :: Python :: 3 :: Only", - "Programming Language :: Python :: 3.6", - "Programming Language :: Python :: 3.7", "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", "Topic :: Scientific/Engineering", "Topic :: Software Development", "Topic :: Utilities", ] dependencies = [ "numpy>=1.14.5", - "pentapy>=1.1.0", + "pentapy>=1.1.0,<2", "scipy>=1.1.0", ] [project.optional-dependencies] doc = [ - "m2r2>=0.2.8", - "matplotlib>=3", + "myst-parser", + "matplotlib>=3.7", "numpydoc>=1.1", - "sphinx>=4", + "sphinx>=7", "sphinx-gallery>=0.8", - "sphinx-rtd-theme>=1,<1.1", + "sphinx-rtd-theme>=3", ] test = ["pytest-cov>=3"] check = [ - "black>=23,<24", - "isort[colors]<6", - "pylint<3", + "black>=25", + "isort[colors]", + "pylint", ] [project.urls] @@ -67,24 +69,39 @@ Tracker = "https://github.com/GeoStat-Framework/anaflow/issues" Changelog = "https://github.com/GeoStat-Framework/anaflow/blob/main/CHANGELOG.md" Conda-Forge = "https://anaconda.org/conda-forge/anaflow" -[tool.setuptools] -license-files = ["LICENSE"] +[tool.hatch.version] +source = "vcs" +fallback_version = "0.0.0.dev0" -[tool.setuptools_scm] -write_to = "src/anaflow/_version.py" -write_to_template = "__version__ = '{version}'" +[tool.hatch.version.raw-options] local_scheme = "no-local-version" -fallback_version = "0.0.0.dev0" + +[tool.hatch.build.hooks.vcs] +version-file = "src/anaflow/_version.py" +template = "__version__ = '{version}'" + +[tool.hatch.build.targets.sdist] +include = [ + "/src", + "/tests", +] + +[tool.hatch.build.targets.wheel] +packages = ["src/anaflow"] [tool.isort] profile = "black" multi_line_output = 3 -line_length = 79 [tool.black] -exclude = "_version.py" -line-length = 79 -target-version = ["py37"] +target-version = [ + "py38", + "py39", + "py310", + "py311", + "py312", + "py313", +] [tool.coverage] [tool.coverage.run] @@ -120,9 +137,10 @@ target-version = ["py37"] output-format = "colorized" [tool.pylint.design] - max-args = 20 + max-args = 25 max-locals = 50 max-branches = 30 max-statements = 80 max-attributes = 25 max-public-methods = 75 + max-positional-arguments=25 diff --git a/src/anaflow/__init__.py b/src/anaflow/__init__.py index 3f22418..f209eb7 100644 --- a/src/anaflow/__init__.py +++ b/src/anaflow/__init__.py @@ -39,10 +39,14 @@ ext_thiem_3d ext_thiem_tpl ext_thiem_tpl_3d + ext_thiem_int + ext_thiem_int_3d ext_theis_2d ext_theis_3d ext_theis_tpl ext_thiem_tpl_3d + ext_theis_int + ext_thiem_int_3d neuman2004 neuman2004_steady @@ -78,15 +82,21 @@ specialrange specialrange_cut """ + +from anaflow import flow, tools from anaflow.flow import ( ext_grf, ext_grf_steady, ext_theis_2d, ext_theis_3d, + ext_theis_int, + ext_theis_int_3d, ext_theis_tpl, ext_theis_tpl_3d, ext_thiem_2d, ext_thiem_3d, + ext_thiem_int, + ext_thiem_int_3d, ext_thiem_tpl, ext_thiem_tpl_3d, grf, @@ -95,13 +105,7 @@ theis, thiem, ) -from anaflow.tools import ( - get_lap, - get_lap_inv, - specialrange, - specialrange_cut, - step_f, -) +from anaflow.tools import get_lap, get_lap_inv, specialrange, specialrange_cut, step_f try: from anaflow._version import __version__ @@ -110,6 +114,7 @@ __version__ = "0.0.0.dev0" __all__ = ["__version__"] +__all__ += ["flow", "tools"] __all__ += [ "thiem", "theis", @@ -117,10 +122,14 @@ "ext_thiem_3d", "ext_thiem_tpl", "ext_thiem_tpl_3d", + "ext_thiem_int", + "ext_thiem_int_3d", "ext_theis_2d", "ext_theis_3d", "ext_theis_tpl", "ext_theis_tpl_3d", + "ext_theis_int", + "ext_theis_int_3d", "neuman2004", "neuman2004_steady", "grf", diff --git a/src/anaflow/flow/__init__.py b/src/anaflow/flow/__init__.py index 059bb19..417f3e8 100644 --- a/src/anaflow/flow/__init__.py +++ b/src/anaflow/flow/__init__.py @@ -38,10 +38,14 @@ ext_thiem_3d ext_thiem_tpl ext_thiem_tpl_3d + ext_thiem_int + ext_thiem_int_3d ext_theis_2d ext_theis_3d ext_theis_tpl ext_theis_tpl_3d + ext_theis_int + ext_theis_int_3d neuman2004 neuman2004_steady @@ -56,14 +60,19 @@ ext_grf ext_grf_steady """ + from anaflow.flow.ext_grf_model import ext_grf, ext_grf_steady from anaflow.flow.heterogeneous import ( ext_theis_2d, ext_theis_3d, + ext_theis_int, + ext_theis_int_3d, ext_theis_tpl, ext_theis_tpl_3d, ext_thiem_2d, ext_thiem_3d, + ext_thiem_int, + ext_thiem_int_3d, ext_thiem_tpl, ext_thiem_tpl_3d, neuman2004, @@ -78,10 +87,14 @@ "ext_thiem_3d", "ext_thiem_tpl", "ext_thiem_tpl_3d", + "ext_thiem_int", + "ext_thiem_int_3d", "ext_theis_2d", "ext_theis_3d", "ext_theis_tpl", "ext_theis_tpl_3d", + "ext_theis_int", + "ext_theis_int_3d", "neuman2004", "neuman2004_steady", "grf", diff --git a/src/anaflow/flow/ext_grf_model.py b/src/anaflow/flow/ext_grf_model.py index 207c8e3..cd1e744 100644 --- a/src/anaflow/flow/ext_grf_model.py +++ b/src/anaflow/flow/ext_grf_model.py @@ -9,6 +9,7 @@ ext_grf ext_grf_steady """ + # pylint: disable=C0103 import numpy as np from scipy.integrate import quad as integ @@ -200,9 +201,7 @@ def integrand(val): if np.isclose(dim, 2): res = np.log(r_ref / Input.rad) / con else: - res = ( - (r_ref ** (2 - dim) - Input.rad ** (2 - dim)) / (2 - dim) / con - ) + res = (r_ref ** (2 - dim) - Input.rad ** (2 - dim)) / (2 - dim) / con res = Input.reshape(res) # rescale by pumping rate diff --git a/src/anaflow/flow/heterogeneous.py b/src/anaflow/flow/heterogeneous.py index 0cfd06b..613dfbf 100644 --- a/src/anaflow/flow/heterogeneous.py +++ b/src/anaflow/flow/heterogeneous.py @@ -10,13 +10,18 @@ ext_thiem_3d ext_thiem_tpl ext_thiem_tpl_3d + ext_thiem_int + ext_thiem_int_3d ext_theis_2d ext_theis_3d ext_theis_tpl ext_theis_tpl_3d + ext_theis_int + ext_theis_int_3d neuman2004 neuman2004_steady """ + # pylint: disable=C0103,C0302 import functools as ft @@ -28,6 +33,8 @@ K_CG, T_CG, TPL_CG, + Int_CG, + Int_CG_error, K_CG_error, T_CG_error, TPL_CG_error, @@ -40,10 +47,14 @@ "ext_thiem_3d", "ext_thiem_tpl", "ext_thiem_tpl_3d", + "ext_thiem_int", + "ext_thiem_int_3d", "ext_theis_2d", "ext_theis_3d", "ext_theis_tpl", "ext_theis_tpl_3d", + "ext_theis_int", + "ext_theis_int_3d", "neuman2004", "neuman2004_steady", ] @@ -93,7 +104,7 @@ def ext_thiem_2d( Reference head at the reference-radius `r_ref`. Default: ``0.0`` T_well : :class:`float`, optional Explicit transmissivity value at the well. Default: ``None`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` @@ -124,13 +135,9 @@ def ext_thiem_2d( rad = np.array(rad, dtype=float) # check the input if r_ref <= 0.0: - raise ValueError( - "The upper boundary needs to be greater than the wellradius" - ) + raise ValueError("The upper boundary needs to be greater than the wellradius") if np.any(rad <= 0.0): - raise ValueError( - "The given radii need to be greater than the wellradius" - ) + raise ValueError("The given radii need to be greater than the wellradius") if trans_gmean <= 0.0: raise ValueError("The Transmissivity needs to be positive.") if T_well is not None and T_well <= 0.0: @@ -209,7 +216,7 @@ def ext_thiem_3d( Explicit conductivity value at the well. One can choose between the harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an arbitrary float value. Default: ``"KH"`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` @@ -313,6 +320,7 @@ def ext_theis_2d( struc_grid=True, far_err=0.01, parts=30, + fix_T_well=False, lap_kwargs=None, ): """ @@ -349,7 +357,7 @@ def ext_theis_2d( Default: ``0.0`` T_well : :class:`float`, optional Explicit transmissivity value at the well. Harmonic mean by default. - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` far_err : :class:`float`, optional @@ -364,6 +372,8 @@ def ext_theis_2d( Since the solution is calculated by setting the transmissivity to local constant values, one needs to specify the number of partitions of the transmissivity. Default: ``30`` + fix_T_well : :class:`bool`, optional + Whether to fix the T_well value in the laplace solver. Default: False lap_kwargs : :class:`dict` or :any:`None` optional Dictionary for :any:`get_lap_inv` containing `method` and `method_dict`. The default is equivalent to @@ -427,7 +437,10 @@ def ext_theis_2d( T_well=T_well, prop=prop, ) - T_well = T_CG(r_well, trans_gmean, var, len_scale, T_well, prop) + if fix_T_well: + T_well = T_CG(r_well, trans_gmean, var, len_scale, T_well, prop) + else: + T_well = None return ext_grf( time=time, rad=rad, @@ -462,6 +475,7 @@ def ext_theis_3d( far_err=0.01, struc_grid=True, parts=30, + fix_K_well=False, lap_kwargs=None, ): """ @@ -506,7 +520,7 @@ def ext_theis_3d( Explicit conductivity value at the well. One can choose between the harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an arbitrary float value. Default: ``"KH"`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` far_err : :class:`float`, optional @@ -521,6 +535,8 @@ def ext_theis_3d( Since the solution is calculated by setting the transmissivity to local constant values, one needs to specify the number of partitions of the transmissivity. Default: ``30`` + fix_K_well : :class:`bool`, optional + Whether to fix the K_well value in the laplace solver. Default: False lap_kwargs : :class:`dict` or :any:`None` optional Dictionary for :any:`get_lap_inv` containing `method` and `method_dict`. The default is equivalent to @@ -569,13 +585,9 @@ def ext_theis_3d( if parts <= 1: raise ValueError("The numbor of partitions needs to be at least 2") if not 0.0 < far_err < 1.0: - raise ValueError( - "The relative error of Conductivity needs to be within (0,1)" - ) + raise ValueError("The relative error of Conductivity needs to be within (0,1)") # genearte rlast from a given relativ-error to farfield-conductivity - r_last = K_CG_error( - far_err, cond_gmean, var, len_scale, anis, K_well, prop - ) + r_last = K_CG_error(far_err, cond_gmean, var, len_scale, anis, K_well, prop) # generate the partition points if r_last > r_well: R_part = specialrange_cut(r_well, r_bound, parts + 1, r_last) @@ -592,7 +604,10 @@ def ext_theis_3d( K_well=K_well, prop=prop, ) - K_well = K_CG(r_well, cond_gmean, var, len_scale, anis, K_well, prop) + if fix_K_well: + K_well = K_CG(r_well, cond_gmean, var, len_scale, anis, K_well, prop) + else: + K_well = None return ext_grf( time=time, rad=rad, @@ -634,6 +649,7 @@ def ext_theis_tpl( far_err=0.01, struc_grid=True, parts=30, + fix_K_well=False, lap_kwargs=None, ): """ @@ -694,7 +710,7 @@ def ext_theis_tpl( Explicit conductivity value at the well. One can choose between the harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an arbitrary float value. Default: ``"KH"`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` far_err : :class:`float`, optional @@ -709,6 +725,8 @@ def ext_theis_tpl( Since the solution is calculated by setting the transmissivity to local constant values, one needs to specify the number of partitions of the transmissivity. Default: ``30`` + fix_K_well : :class:`bool`, optional + Whether to fix the K_well value in the laplace solver. Default: False lap_kwargs : :class:`dict` or :any:`None` optional Dictionary for :any:`get_lap_inv` containing `method` and `method_dict`. The default is equivalent to @@ -753,9 +771,7 @@ def ext_theis_tpl( if parts <= 1: raise ValueError("The numbor of partitions needs to be at least 2") if not 0.0 < far_err < 1.0: - raise ValueError( - "The relative error of Conductivity needs to be within (0,1)" - ) + raise ValueError("The relative error of Conductivity needs to be within (0,1)") # genearte rlast from a given relativ-error to farfield-conductivity r_last = TPL_CG_error( far_err, cond_gmean, len_scale, hurst, var, c, 1, dim, K_well, prop @@ -780,9 +796,12 @@ def ext_theis_tpl( K_well=K_well, prop=prop, ) - K_well = TPL_CG( - r_well, cond_gmean, len_scale, hurst, var, c, 1, dim, K_well, prop - ) + if fix_K_well: + K_well = TPL_CG( + r_well, cond_gmean, len_scale, hurst, var, c, 1, dim, K_well, prop + ) + else: + K_well = None return ext_grf( time=time, rad=rad, @@ -819,6 +838,7 @@ def ext_theis_tpl_3d( far_err=0.01, struc_grid=True, parts=30, + fix_K_well=False, lap_kwargs=None, ): """ @@ -873,7 +893,7 @@ def ext_theis_tpl_3d( Explicit conductivity value at the well. One can choose between the harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an arbitrary float value. Default: ``"KH"`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` far_err : :class:`float`, optional @@ -888,6 +908,8 @@ def ext_theis_tpl_3d( Since the solution is calculated by setting the transmissivity to local constant values, one needs to specify the number of partitions of the transmissivity. Default: ``30`` + fix_K_well : :class:`bool`, optional + Whether to fix the K_well value in the laplace solver. Default: False lap_kwargs : :class:`dict` or :any:`None` optional Dictionary for :any:`get_lap_inv` containing `method` and `method_dict`. The default is equivalent to @@ -932,9 +954,7 @@ def ext_theis_tpl_3d( if parts <= 1: raise ValueError("The numbor of partitions needs to be at least 2") if not 0.0 < far_err < 1.0: - raise ValueError( - "The relative error of Conductivity needs to be within (0,1)" - ) + raise ValueError("The relative error of Conductivity needs to be within (0,1)") # genearte rlast from a given relativ-error to farfield-conductivity r_last = TPL_CG_error( far_err, cond_gmean, len_scale, hurst, var, c, anis, 3, K_well, prop @@ -959,9 +979,12 @@ def ext_theis_tpl_3d( K_well=K_well, prop=prop, ) - K_well = TPL_CG( - r_well, cond_gmean, len_scale, hurst, var, c, anis, 3, K_well, prop - ) + if fix_K_well: + K_well = TPL_CG( + r_well, cond_gmean, len_scale, hurst, var, c, anis, 3, K_well, prop + ) + else: + K_well = None return ext_grf( time=time, rad=rad, @@ -1044,7 +1067,7 @@ def ext_thiem_tpl( Explicit conductivity value at the well. One can choose between the harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an arbitrary float value. Default: ``"KH"`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` @@ -1167,7 +1190,7 @@ def ext_thiem_tpl_3d( Explicit conductivity value at the well. One can choose between the harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an arbitrary float value. Default: ``"KH"`` - prop: :class:`float`, optional + prop : :class:`float`, optional Proportionality factor used within the upscaling procedure. Default: ``1.6`` far_err : :class:`float`, optional @@ -1241,47 +1264,76 @@ def ext_thiem_tpl_3d( ############################################################################### -# solution for apparent transmissivity from Neuman 2004 +# extended Thiem/Theis solutions for the Integral model ############################################################################### -def neuman2004( +def ext_theis_int( time, rad, storage, - trans_gmean, + cond_gmean, var, len_scale, + roughness=1.0, + dim=2.0, + lat_ext=1.0, rate=-1e-4, r_well=0.0, r_bound=np.inf, h_bound=0.0, + K_well="KH", + prop=1.6, + far_err=0.01, struc_grid=True, parts=30, + fix_K_well=False, lap_kwargs=None, ): """ - The transient solution for the apparent transmissivity from [Neuman2004]. + The extended Theis solution for the integral variogram model. - This solution is build on the apparent transmissivity from Neuman 2004, - which represents a mean drawdown in an ensemble of pumping tests in - heterogeneous transmissivity fields following an exponential covariance. - Presented in [Neuman2004]_. + The extended Theis solution for transient flow under + a pumping condition in a confined aquifer. + The type curve is describing the effective drawdown + in a d-dimensional statistical framework, + where the conductivity distribution is + following a log-normal distribution with an integral + correlation incorporating a roughness parameter. + + The roughness parameter controls the field roughness of + log-transmissivity. It's limits are a pure nugget model for + r -> 0 and the gaussian model for r -> infinity. Parameters ---------- time : :class:`numpy.ndarray` - Array with all time-points where the function should be evaluated. + 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. + Array with all radii where the function should be evaluated storage : :class:`float` Storage of the aquifer. - trans_gmean : :class:`float` - Geometric-mean transmissivity. + cond_gmean : :class:`float` + Geometric-mean conductivity. + You can also treat this as transmissivity by leaving 'lat_ext=1'. var : :class:`float` - Variance of log-transmissivity. + Variance of the log-conductivity. len_scale : :class:`float` - Correlation-length of log-transmissivity. + Corralation-length of log-conductivity. + roughness: :class:`float`, optional + Roughness of the model. Should be positive. + Default: ``1`` + dim: :class:`float`, optional + Dimension of space. + Default: ``2.0`` + lat_ext : :class:`float`, optional + Lateral extend of the aquifer: + + * sqare-root of cross-section in 1D + * thickness in 2D + * meaningless in 3D + + Default: ``1.0`` rate : :class:`float`, optional Pumpingrate at the well. Default: -1e-4 r_well : :class:`float`, optional @@ -1291,6 +1343,16 @@ def neuman2004( h_bound : :class:`float`, optional Reference head at the outer boundary as well as initial condition. Default: ``0.0`` + K_well : :class:`float`, optional + Explicit conductivity value at the well. One can choose between the + harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an + arbitrary float value. Default: ``"KH"`` + prop : :class:`float`, optional + Proportionality factor used within the upscaling procedure. + Default: ``1.6`` + far_err : :class:`float`, optional + Relative error for the farfield transmissivity for calculating the + cutoff-point of the solution. Default: ``0.01`` 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 @@ -1300,6 +1362,8 @@ def neuman2004( Since the solution is calculated by setting the transmissivity to local constant values, one needs to specify the number of partitions of the transmissivity. Default: ``30`` + fix_K_well : :class:`bool`, optional + Whether to fix the K_well value in the laplace solver. Default: False lap_kwargs : :class:`dict` or :any:`None` optional Dictionary for :any:`get_lap_inv` containing `method` and `method_dict`. The default is equivalent to @@ -1311,11 +1375,10 @@ def neuman2004( 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 + Notes + ----- + If you want to use cartesian coordiantes, just use the formula + ``r = sqrt(x**2 + y**2)`` """ # check the input if r_well < 0.0: @@ -1323,50 +1386,631 @@ def neuman2004( if r_bound <= r_well: raise ValueError("The upper boundary needs to be > well radius") if storage <= 0.0: - raise ValueError("The Storage needs to be positive.") - if trans_gmean <= 0.0: - raise ValueError("The Transmissivity needs to be positive.") - if var < 0.0: - raise ValueError("The variance needs to be positive.") + raise ValueError("The storage needs to be positive.") + if cond_gmean <= 0.0: + raise ValueError("The gmean conductivity needs to be positive.") if len_scale <= 0.0: raise ValueError("The correlationlength needs to be positive.") + if roughness <= 0.0: + raise ValueError("Roughness needs to be positive.") + if var < 0.0: + raise ValueError("The variance needs to be positive.") + if K_well != "KA" and K_well != "KH" and not isinstance(K_well, float): + raise ValueError( + "The well-conductivity should be given as float or 'KA' resp 'KH'" + ) + if isinstance(K_well, float) and K_well <= 0.0: + raise ValueError("The well-conductivity needs to be positive.") + if prop <= 0.0: + raise ValueError("The proportionality factor needs to be positive.") 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 - r_last = 2 * len_scale + if not 0.0 < far_err < 1.0: + raise ValueError("The relative error of Conductivity needs to be within (0,1)") + # genearte rlast from a given relativ-error to farfield-conductivity + r_last = Int_CG_error( + far_err, cond_gmean, var, len_scale, roughness, 1, dim, K_well, prop + ) # generate the partition points if r_last > r_well: R_part = specialrange_cut(r_well, r_bound, parts + 1, r_last) else: R_part = np.array([r_well, r_bound]) - # calculate the harmonic mean transmissivity values within each partition - T_part = annular_hmean( - neuman2004_trans, + # calculate the harmonic mean conductivity values within each partition + K_part = annular_hmean( + Int_CG, R_part, - trans_gmean=trans_gmean, + ann_dim=dim, + cond_gmean=cond_gmean, var=var, len_scale=len_scale, + roughness=roughness, + anis=1, + dim=dim, + K_well=K_well, + prop=prop, ) - T_well = neuman2004_trans(r_well, trans_gmean, var, len_scale) + if fix_K_well: + K_well = Int_CG( + r_well, cond_gmean, var, len_scale, roughness, 1, dim, K_well, prop + ) + else: + K_well = None return ext_grf( time=time, rad=rad, - S_part=np.full_like(T_part, storage), - K_part=T_part, + S_part=np.full_like(K_part, storage), + K_part=K_part, R_part=R_part, - dim=2, - lat_ext=1, + dim=dim, + lat_ext=lat_ext, rate=rate, h_bound=h_bound, - K_well=T_well, + K_well=K_well, struc_grid=struc_grid, lap_kwargs=lap_kwargs, ) -def neuman2004_steady( - rad, r_ref, trans_gmean, var, len_scale, rate=-1e-4, h_ref=0.0 +def ext_theis_int_3d( + time, + rad, + storage, + cond_gmean, + var, + len_scale, + roughness=1.0, + anis=1, + lat_ext=1.0, + rate=-1e-4, + r_well=0.0, + r_bound=np.inf, + h_bound=0.0, + K_well="KH", + prop=1.6, + far_err=0.01, + struc_grid=True, + parts=30, + fix_K_well=False, + lap_kwargs=None, ): + """ + The extended Theis solution for the integral variogram in 3D. + + The extended Theis solution for transient flow under + a pumping condition in a confined aquifer with anisotropy in 3D. + The type curve is describing the effective drawdown + in a 3-dimensional statistical framework, + where the conductivity distribution is + following a log-normal distribution with an Integral + correlation function incorporating a roughness parameter. + + The roughness parameter controls the field roughness of + log-transmissivity. It's limits are a pure nugget model for + r -> 0 and the gaussian model for r -> infinity. + + 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` + Storage of the aquifer. + cond_gmean : :class:`float` + Geometric-mean conductivity. + var : :class:`float` + Variance of the log-conductivity. + len_scale : :class:`float` + Corralation-length of log-conductivity. + roughness: :class:`float`, optional + Roughness of the model. Should be positive. + Default: ``1`` + anis : :class:`float`, optional + Anisotropy-ratio of the vertical and horizontal corralation-lengths. + Default: 1.0 + lat_ext : :class:`float`, optional + Lateral extend of the aquifer (thickness). + Default: ``1.0`` + rate : :class:`float`, optional + Pumpingrate at the well. Default: -1e-4 + r_well : :class:`float`, optional + 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`` + K_well : :class:`float`, optional + Explicit conductivity value at the well. One can choose between the + harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an + arbitrary float value. Default: ``"KH"`` + prop : :class:`float`, optional + Proportionality factor used within the upscaling procedure. + Default: ``1.6`` + far_err : :class:`float`, optional + Relative error for the farfield transmissivity for calculating the + cutoff-point of the solution. Default: ``0.01`` + 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`` + parts : :class:`int`, optional + Since the solution is calculated by setting the transmissivity to local + constant values, one needs to specify the number of partitions of the + transmissivity. Default: ``30`` + fix_K_well : :class:`bool`, optional + Whether to fix the K_well value in the laplace solver. Default: False + 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. + + Notes + ----- + If you want to use cartesian coordiantes, just use the formula + ``r = sqrt(x**2 + y**2)`` + """ + # 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 > well radius") + if storage <= 0.0: + raise ValueError("The storage needs to be positive.") + if cond_gmean <= 0.0: + raise ValueError("The gmean conductivity needs to be positive.") + if len_scale <= 0.0: + raise ValueError("The correlationlength needs to be positive.") + if roughness <= 0.0: + raise ValueError("Roughness needs to be positive.") + if var < 0.0: + raise ValueError("The variance needs to be positive.") + if K_well != "KA" and K_well != "KH" and not isinstance(K_well, float): + raise ValueError( + "The well-conductivity should be given as float or 'KA' resp 'KH'" + ) + if isinstance(K_well, float) and K_well <= 0.0: + raise ValueError("The well-conductivity needs to be positive.") + if prop <= 0.0: + raise ValueError("The proportionality factor needs to be positive.") + if parts <= 1: + raise ValueError("The numbor of partitions needs to be at least 2") + if not 0.0 < far_err < 1.0: + raise ValueError("The relative error of Conductivity needs to be within (0,1)") + # genearte rlast from a given relativ-error to farfield-conductivity + r_last = Int_CG_error( + far_err, cond_gmean, var, len_scale, roughness, anis, 3, K_well, prop + ) + # generate the partition points + if r_last > r_well: + R_part = specialrange_cut(r_well, r_bound, parts + 1, r_last) + else: + R_part = np.array([r_well, r_bound]) + # calculate the harmonic mean conductivity values within each partition + K_part = annular_hmean( + Int_CG, + R_part, + ann_dim=2, + cond_gmean=cond_gmean, + var=var, + len_scale=len_scale, + roughness=roughness, + anis=anis, + dim=3, + K_well=K_well, + prop=prop, + ) + if fix_K_well: + K_well = Int_CG( + r_well, cond_gmean, var, len_scale, roughness, anis, 3, K_well, prop + ) + else: + K_well = None + return ext_grf( + time=time, + rad=rad, + S_part=np.full_like(K_part, storage), + K_part=K_part, + R_part=R_part, + dim=2, + lat_ext=lat_ext, + rate=rate, + h_bound=h_bound, + K_well=K_well, + struc_grid=struc_grid, + lap_kwargs=lap_kwargs, + ) + + +def ext_thiem_int( + rad, + r_ref, + cond_gmean, + var, + len_scale, + roughness=1.0, + dim=2.0, + lat_ext=1.0, + rate=-1e-4, + h_ref=0.0, + K_well="KH", + prop=1.6, +): + """ + The extended Thiem solution for the intgeral variogram. + + The extended Theis solution for steady flow under + a pumping condition in a confined aquifer. + The type curve is describing the effective drawdown + in a d-dimensional statistical framework, + where the conductivity distribution is + following a log-normal distribution with an integral + correlation function incorporating a roughness parameter. + + The roughness parameter controls the field roughness of + log-transmissivity. It's limits are a pure nugget model for + r -> 0 and the gaussian model for r -> infinity. + + Parameters + ---------- + rad : :class:`numpy.ndarray` + Array with all radii where the function should be evaluated + r_ref : :class:`float` + Reference radius with known head (see `h_ref`) + cond_gmean : :class:`float` + Geometric-mean conductivity. + You can also treat this as transmissivity by leaving 'lat_ext=1'. + var : :class:`float` + Variance of the log-conductivity. + len_scale : :class:`float` + Corralation-length of log-conductivity. + roughness: :class:`float`, optional + Roughness of the model. Should be positive. + Default: ``1`` + dim: :class:`float`, optional + Dimension of space. + Default: ``2.0`` + lat_ext : :class:`float`, optional + Lateral extend of the aquifer: + + * sqare-root of cross-section in 1D + * thickness in 2D + * meaningless in 3D + + Default: ``1.0`` + rate : :class:`float`, optional + Pumpingrate at the well. Default: -1e-4 + h_ref : :class:`float`, optional + Reference head at the reference-radius `r_ref`. Default: ``0.0`` + K_well : :class:`float`, optional + Explicit conductivity value at the well. One can choose between the + harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an + arbitrary float value. Default: ``"KH"`` + prop : :class:`float`, optional + Proportionality factor used within the upscaling procedure. + Default: ``1.6`` + + Returns + ------- + head : :class:`numpy.ndarray` + Array with all heads at the given radii and time-points. + + Notes + ----- + If you want to use cartesian coordiantes, just use the formula + ``r = sqrt(x**2 + y**2)`` + """ + # check the input + if cond_gmean <= 0.0: + raise ValueError("The gmean conductivity needs to be positive.") + if len_scale <= 0.0: + raise ValueError("The correlationlength needs to be positive.") + if roughness <= 0.0: + raise ValueError("Roughness needs to be positive.") + if var < 0.0: + raise ValueError("The variance needs to be positive.") + if K_well != "KA" and K_well != "KH" and not isinstance(K_well, float): + raise ValueError( + "The well-conductivity should be given as float or 'KA' resp 'KH'" + ) + if isinstance(K_well, float) and K_well <= 0.0: + raise ValueError("The well-conductivity needs to be positive.") + if prop <= 0.0: + raise ValueError("The proportionality factor needs to be positive.") + cond = ft.partial( + Int_CG, + cond_gmean=cond_gmean, + len_scale=len_scale, + var=var, + roughness=roughness, + anis=1, + dim=dim, + K_well=K_well, + prop=prop, + ) + return ext_grf_steady( + rad=rad, + r_ref=r_ref, + conductivity=cond, + dim=dim, + lat_ext=lat_ext, + rate=rate, + h_ref=h_ref, + ) + + +def ext_thiem_int_3d( + rad, + r_ref, + cond_gmean, + var, + len_scale, + roughness=1.0, + anis=1, + lat_ext=1.0, + rate=-1e-4, + h_ref=0.0, + K_well="KH", + prop=1.6, +): + """ + The extended Theis solution for the integral variogram in 3D. + + The extended Theis solution for transient flow under + a pumping condition in a confined aquifer with anisotropy in 3D. + The type curve is describing the effective drawdown + in a 3-dimensional statistical framework, + where the conductivity distribution is + following a log-normal distribution with an integral + correlation function incorporating a roughness parameter. + + The roughness parameter controls the field roughness of + log-transmissivity. It's limits are a pure nugget model for + r -> 0 and the gaussian model for r -> infinity. + + 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` + Storage of the aquifer. + cond_gmean : :class:`float` + Geometric-mean conductivity. + var : :class:`float` + Variance of the log-conductivity. + len_scale : :class:`float` + Corralation-length of log-conductivity. + roughness: :class:`float`, optional + Roughness of the model. Should be positive. + Default: ``1`` + anis : :class:`float`, optional + Anisotropy-ratio of the vertical and horizontal corralation-lengths. + Default: 1.0 + lat_ext : :class:`float`, optional + Lateral extend of the aquifer (thickness). + Default: ``1.0`` + rate : :class:`float`, optional + Pumpingrate at the well. Default: -1e-4 + r_well : :class:`float`, optional + 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`` + K_well : :class:`float`, optional + Explicit conductivity value at the well. One can choose between the + harmonic mean (``"KH"``), the arithmetic mean (``"KA"``) or an + arbitrary float value. Default: ``"KH"`` + prop : :class:`float`, optional + Proportionality factor used within the upscaling procedure. + Default: ``1.6`` + far_err : :class:`float`, optional + Relative error for the farfield transmissivity for calculating the + cutoff-point of the solution. Default: ``0.01`` + 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`` + parts : :class:`int`, optional + Since the solution is calculated by setting the transmissivity 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. + + Notes + ----- + If you want to use cartesian coordiantes, just use the formula + ``r = sqrt(x**2 + y**2)`` + """ + # check the input + if cond_gmean <= 0.0: + raise ValueError("The gmean conductivity needs to be positive.") + if len_scale <= 0.0: + raise ValueError("The correlationlength needs to be positive.") + if roughness <= 0.0: + raise ValueError("Roughness needs to be positive.") + if var < 0.0: + raise ValueError("The variance needs to be positive.") + if K_well != "KA" and K_well != "KH" and not isinstance(K_well, float): + raise ValueError( + "The well-conductivity should be given as float or 'KA' resp 'KH'" + ) + if isinstance(K_well, float) and K_well <= 0.0: + raise ValueError("The well-conductivity needs to be positive.") + if prop <= 0.0: + raise ValueError("The proportionality factor needs to be positive.") + cond = ft.partial( + Int_CG, + cond_gmean=cond_gmean, + var=var, + len_scale=len_scale, + roughness=roughness, + anis=anis, + dim=3, + K_well=K_well, + prop=prop, + ) + return ext_grf_steady( + rad=rad, + r_ref=r_ref, + conductivity=cond, + dim=2, + lat_ext=lat_ext, + rate=rate, + h_ref=h_ref, + ) + + +############################################################################### +# solution for apparent transmissivity from Neuman 2004 +############################################################################### + + +def neuman2004( + time, + rad, + storage, + trans_gmean, + var, + len_scale, + rate=-1e-4, + r_well=0.0, + r_bound=np.inf, + h_bound=0.0, + struc_grid=True, + parts=30, + fix_T_well=False, + lap_kwargs=None, +): + """ + The transient solution for the apparent transmissivity from [Neuman2004]. + + This solution is build on the apparent transmissivity from Neuman 2004, + which represents a mean drawdown in an ensemble of pumping tests in + heterogeneous transmissivity fields following an exponential covariance. + Presented in [Neuman2004]_. + + 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` + Storage 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 + r_well : :class:`float`, optional + 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`` + 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`` + parts : :class:`int`, optional + Since the solution is calculated by setting the transmissivity 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 + """ + # 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 > well radius") + if storage <= 0.0: + raise ValueError("The Storage needs to be positive.") + if trans_gmean <= 0.0: + raise ValueError("The Transmissivity needs to be positive.") + if var < 0.0: + raise ValueError("The variance needs to be positive.") + if len_scale <= 0.0: + raise ValueError("The correlationlength needs to be positive.") + 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 + r_last = 2 * len_scale + # generate the partition points + if r_last > r_well: + R_part = specialrange_cut(r_well, r_bound, parts + 1, r_last) + else: + R_part = np.array([r_well, r_bound]) + # calculate the harmonic mean transmissivity values within each partition + T_part = annular_hmean( + neuman2004_trans, + R_part, + trans_gmean=trans_gmean, + var=var, + len_scale=len_scale, + ) + if fix_T_well: + T_well = neuman2004_trans(r_well, trans_gmean, var, len_scale) + else: + T_well = None + return ext_grf( + time=time, + rad=rad, + S_part=np.full_like(T_part, storage), + K_part=T_part, + R_part=R_part, + dim=2, + lat_ext=1, + rate=rate, + h_bound=h_bound, + K_well=T_well, + struc_grid=struc_grid, + lap_kwargs=lap_kwargs, + ) + + +def neuman2004_steady(rad, r_ref, trans_gmean, var, len_scale, rate=-1e-4, h_ref=0.0): """ The steady solution for the apparent transmissivity from [Neuman2004]. diff --git a/src/anaflow/flow/homogeneous.py b/src/anaflow/flow/homogeneous.py index 2fa0865..f1f878c 100644 --- a/src/anaflow/flow/homogeneous.py +++ b/src/anaflow/flow/homogeneous.py @@ -10,6 +10,7 @@ theis grf """ + # pylint: disable=C0103 import numpy as np diff --git a/src/anaflow/flow/laplace.py b/src/anaflow/flow/laplace.py index a7256ce..a47b917 100644 --- a/src/anaflow/flow/laplace.py +++ b/src/anaflow/flow/laplace.py @@ -10,6 +10,7 @@ grf_laplace """ + # pylint: disable=C0103,R0915 import warnings @@ -200,9 +201,7 @@ def grf_laplace( # calculate the head for ri, re in enumerate(rad): if re < R_part[-1]: - res[si, ri] = re**nu * ( - As * kv(nu, Cs * re) + Bs * iv(nu, Cs * re) - ) + res[si, ri] = re**nu * (As * kv(nu, Cs * re) + Bs * iv(nu, Cs * re)) # if there is more than one partition, create an equation system else: diff --git a/src/anaflow/tools/__init__.py b/src/anaflow/tools/__init__.py index 6d9671d..8a3a0dd 100644 --- a/src/anaflow/tools/__init__.py +++ b/src/anaflow/tools/__init__.py @@ -42,6 +42,7 @@ T_CG K_CG TPL_CG + Int_CG Special ~~~~~~~ @@ -68,7 +69,8 @@ get_lap get_lap_inv """ -from anaflow.tools.coarse_graining import K_CG, T_CG, TPL_CG + +from anaflow.tools.coarse_graining import K_CG, T_CG, TPL_CG, Int_CG from anaflow.tools.laplace import get_lap, get_lap_inv from anaflow.tools.mean import ( annular_amean, @@ -101,4 +103,5 @@ "T_CG", "K_CG", "TPL_CG", + "Int_CG", ] diff --git a/src/anaflow/tools/coarse_graining.py b/src/anaflow/tools/coarse_graining.py index a2f28ac..2529700 100644 --- a/src/anaflow/tools/coarse_graining.py +++ b/src/anaflow/tools/coarse_graining.py @@ -16,7 +16,10 @@ K_CG_error TPL_CG TPL_CG_error + Int_CG + Int_CG_error """ + # pylint: disable=C0103 import numpy as np from scipy.optimize import root @@ -32,6 +35,8 @@ "K_CG_error", "TPL_CG", "TPL_CG_error", + "Int_CG", + "Int_CG_error", ] @@ -225,9 +230,7 @@ def K_CG(rad, cond_gmean, var, len_scale, anis, K_well="KH", prop=1.6): chi = np.log(K_well / K_efu) return K_efu * np.exp( - chi - / np.sqrt(1.0 + (prop * rad / (len_scale * anis ** (1.0 / 3.0))) ** 2) - ** 3 + chi / np.sqrt(1.0 + (prop * rad / (len_scale * anis ** (1.0 / 3.0))) ** 2) ** 3 ) @@ -334,9 +337,7 @@ def K_CG_error(err, cond_gmean, var, len_scale, anis, K_well="KH", prop=1.6): if chi > 0.0: if chi / np.log(1.0 + err) >= 1.0: - return coef * np.sqrt( - (chi / np.log(1.0 + err)) ** (2.0 / 3.0) - 1.0 - ) + return coef * np.sqrt((chi / np.log(1.0 + err)) ** (2.0 / 3.0) - 1.0) # standard value if the error is less then the variation return 0 @@ -522,6 +523,171 @@ def curve(x): return root(curve, 2 * len_scale)["x"][0] +def Int_CG( + rad, + cond_gmean, + var, + len_scale, + roughness=1.0, + anis=1, + dim=2.0, + K_well="KH", + prop=1.6, +): + """ + The integral coarse-graining conductivity. + + The roughness parameter controls the field roughness of + log-transmissivity. It's limits are a pure nugget model for + r -> 0 and the gaussian model for r -> infinity. + + Parameters + ---------- + rad : :class:`numpy.ndarray` + Array with all radii where the function should be evaluated + cond_gmean : :class:`float` + Geometric-mean conductivity. + var : :class:`float` or :any:`None`, optional + Variance of log-conductivity. + len_scale : :class:`float` + Correlation-length of log-conductivity. + roughness: :class:`float`, optional + Roughness of the model. Should be positive. + Default: 1.0 + anis : :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`` + K_well : :class:`str` or :class:`float`, optional + Explicit conductivity value at the well. One can choose between the + harmonic mean (``"KH"``), + the arithmetic mean (``"KA"``) or an arbitrary float + value. Default: ``"KH"`` + prop: :class:`float`, optional + Proportionality factor used within the upscaling procedure. + Default: ``1.6`` + + Returns + ------- + Int_CG : :class:`numpy.ndarray` + Array containing the effective conductivity values. + """ + # handle special case in 3D with anisotropy + anis = 1 if not np.isclose(dim, 3) else anis + ani = aniso(anis) if np.isclose(dim, 3) else 1 / dim + K_efu = cond_gmean * np.exp(var * (0.5 - ani)) + if K_well == "KH": + chi = var * (ani - 1) + elif K_well == "KA": + chi = var * ani + else: + chi = np.log(K_well / K_efu) + + return K_efu * np.exp( + (chi * roughness / (dim + roughness)) + * tpl_hyp(rad, dim, roughness / 2, len_scale * anis ** (1 / 3), prop) + ) + + +def Int_CG_error( + err, + cond_gmean, + var, + len_scale, + roughness=1.0, + anis=1, + dim=2.0, + K_well="KH", + prop=1.6, +): + """ + Calculating the radial-point for given error. + + Calculating the radial-point where the relative error of the farfield + value is less than the given tollerance. + See: :func:`Int_CG` + + Parameters + ---------- + err : :class:`float` + Given relative error for the farfield conductivity + cond_gmean : :class:`float` + Geometric-mean conductivity + var : :class:`float` or :any:`None`, optional + Variance of log-conductivity. + len_scale : :class:`float` + Correlation-length of log-conductivity. + roughness: :class:`float`, optional + Roughness of the model. Should be positive. + Default: 1.0 + anis : :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`` + K_well : :class:`str` or :class:`float`, optional + Explicit conductivity value at the well. One can choose between the + harmonic mean (``"KH"``), + the arithmetic mean (``"KA"``) or an arbitrary float + value. Default: ``"KH"`` + prop: :class:`float`, optional + Proportionality factor used within the upscaling procedure. + Default: ``1.6`` + + Returns + ------- + rad : :class:`float` + Radial point, where the relative error is less than the given one. + """ + # handle special case in 3D with anisotropy + anis = 1 if not np.isclose(dim, 3) else anis + ani = aniso(anis) if np.isclose(dim, 3) else 1 / dim + K_efu = cond_gmean * np.exp(var * (0.5 - ani)) + if K_well == "KH": + chi = var * (ani - 1) + elif K_well == "KA": + chi = var * ani + else: + chi = np.log(K_well / K_efu) + Kw = np.exp(chi + np.log(K_efu)) + + # define a curve, that has its root at the wanted percentile + if chi > 0: + per = (1 + err) * K_efu + if not per < Kw: + return 0 + elif chi < 0: + per = (1 - err) * K_efu + if not per > Kw: + return 0 + else: + return 0 + + def curve(x): + """Curve for fitting.""" + return ( + Int_CG( + x, + cond_gmean=cond_gmean, + var=var, + len_scale=len_scale, + roughness=roughness, + anis=anis, + dim=dim, + K_well=K_well, + prop=prop, + ) + - per + ) + + return root(curve, 2 * len_scale)["x"][0] + + if __name__ == "__main__": import doctest diff --git a/src/anaflow/tools/laplace.py b/src/anaflow/tools/laplace.py index 467fd78..73d597e 100644 --- a/src/anaflow/tools/laplace.py +++ b/src/anaflow/tools/laplace.py @@ -13,6 +13,7 @@ lap_trans stehfest """ + from math import factorial, floor import numpy as np @@ -131,9 +132,7 @@ def integrand(val): return result -def get_lap_inv( - func, method="stehfest", method_dict=None, arg_dict=None, **kwargs -): +def get_lap_inv(func, method="stehfest", method_dict=None, arg_dict=None, **kwargs): """ Callable Laplace inversion. @@ -295,13 +294,9 @@ def stehfest(func, time, bound=12, arg_dict=None, **kwargs): "The time-values need to be positiv for the stehfest-algorithm" ) if bound <= 1: - raise ValueError( - "The boundary needs to be >1 for the stehfest-algorithm" - ) + raise ValueError("The boundary needs to be >1 for the stehfest-algorithm") if bound % 2 != 0: - raise ValueError( - "The boundary needs to be even for the stehfest-algorithm" - ) + raise ValueError("The boundary needs to be even for the stehfest-algorithm") # get all coefficient factors at once c_fac = c_array(bound) @@ -457,11 +452,7 @@ def _c(i, bound): def _d(k, i, bound): res = ((float(k)) ** (bound / 2 + 1)) * (factorial(2 * k)) - res /= ( - (factorial(bound / 2 - k)) - * (factorial(i - k)) - * (factorial(2 * k - i)) - ) + res /= (factorial(bound / 2 - k)) * (factorial(i - k)) * (factorial(2 * k - i)) res /= factorial(k) ** 2 return res diff --git a/src/anaflow/tools/mean.py b/src/anaflow/tools/mean.py index a06323f..0955f8f 100644 --- a/src/anaflow/tools/mean.py +++ b/src/anaflow/tools/mean.py @@ -14,6 +14,7 @@ annular_hmean annular_pmean """ + # pylint: disable=E1137, C0103 import numpy as np from scipy.integrate import quad as integ @@ -27,9 +28,7 @@ ] -def annular_fmean( - func, val_arr, f_def, f_inv, ann_dim=2, arg_dict=None, **kwargs -): +def annular_fmean(func, val_arr, f_def, f_inv, ann_dim=2, arg_dict=None, **kwargs): r""" Calculating the annular generalized f-mean. @@ -101,9 +100,7 @@ def annular_fmean( if not callable(f_inv): raise ValueError("The inverse f-mean function needs to be callable") if not np.all( - np.isclose( - f_inv(f_def(func(val_arr, **kwargs))), func(val_arr, **kwargs) - ) + np.isclose(f_inv(f_def(func(val_arr, **kwargs))), func(val_arr, **kwargs)) ): raise ValueError("f_def and f_inv need to be inverse to each other") if len(val_arr) < 2: diff --git a/src/anaflow/tools/special.py b/src/anaflow/tools/special.py index 98888fd..a18adb7 100644 --- a/src/anaflow/tools/special.py +++ b/src/anaflow/tools/special.py @@ -20,6 +20,7 @@ tpl_hyp neuman2004_trans """ + # pylint: disable=C0103,R0903 import numpy as np from scipy.special import exp1, expn, gamma, gammaincc, hyp2f1 @@ -77,12 +78,8 @@ def __init__(self, time=0, rad=0, struc_grid=True): self.rad_max = np.max(self.rad) self.time_gz = self.time > 0 - self.time_mat = np.outer( - self.time[self.time_gz], np.ones_like(self.rad) - ) - self.rad_mat = np.outer( - np.ones_like(self.time[self.time_gz]), self.rad - ) + self.time_mat = np.outer(self.time[self.time_gz], np.ones_like(self.rad)) + self.rad_mat = np.outer(np.ones_like(self.time[self.time_gz]), self.rad) if not self.struc_grid and self.rad_shape != self.time_shape: raise ValueError("No struc_grid: shape of time & radius differ") @@ -157,9 +154,7 @@ def specialrange(val_min, val_max, steps, typ="exp"): array([ 1. , 2.53034834, 5.23167968, 10. ]) """ if typ in ["exponential", "exp"]: - rng = np.expm1( - np.linspace(np.log1p(val_min), np.log1p(val_max), steps) - ) + rng = np.expm1(np.linspace(np.log1p(val_min), np.log1p(val_max), steps)) elif typ in ["logarithmic", "log"]: rng = np.log(np.linspace(np.exp(val_min), np.exp(val_max), steps)) elif typ in ["geometric", "geo", "geom"]: @@ -170,9 +165,7 @@ def specialrange(val_min, val_max, steps, typ="exp"): rng = (np.linspace(np.sqrt(val_min), np.sqrt(val_max), steps)) ** 2 elif typ in ["cubic", "cub"]: rng = ( - np.linspace( - np.power(val_min, 1 / 3.0), np.power(val_max, 1 / 3.0), steps - ) + np.linspace(np.power(val_min, 1 / 3.0), np.power(val_max, 1 / 3.0), steps) ) ** 3 elif isinstance(typ, (float, int)): rng = ( @@ -279,12 +272,7 @@ def aniso(e): res = 0.0 else: res = e / (2 * (1.0 - e**2)) - res *= ( - 1.0 - / np.sqrt(1.0 - e**2) - * np.arctan(np.sqrt(1.0 / e**2 - 1.0)) - - e - ) + res *= 1.0 / np.sqrt(1.0 - e**2) * np.arctan(np.sqrt(1.0 / e**2 - 1.0)) - e return res @@ -499,8 +487,8 @@ def inc_gamma(s, x): def tpl_hyp(rad, dim, hurst, corr, prop): """Hyp_2F1 for the TPL CG model.""" - x = 1.0 / (1.0 + (prop * rad / corr) ** 2) - return x ** (dim / 2.0) * hyp2f1(dim / 2.0, 1, dim / 2.0 + 1 + hurst, x) + x, d = 1 / (1 + (prop * rad / corr) ** 2), dim / 2 + return x**d * hyp2f1(d, 1, d + 1 + hurst, x) def neuman2004_trans(rad, trans_gmean, var, len_scale): diff --git a/tests/test_anaflow.py b/tests/test_anaflow.py index afa76d3..7cfcf92 100644 --- a/tests/test_anaflow.py +++ b/tests/test_anaflow.py @@ -63,9 +63,7 @@ def test_homogeneous(self): self.assertTrue(inc(rad_arr, self.delta)) for time_arr in transient.T: self.assertTrue(dec(time_arr, self.delta)) - self.assertAlmostEqual( - np.abs(np.max(steady - transient[-1, :])), 0.0, places=4 - ) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=4) def test_ext_2d(self): steady = af.ext_thiem_2d( @@ -93,9 +91,7 @@ def test_ext_2d(self): self.assertTrue(inc(rad_arr, self.delta)) for time_arr in transient.T: self.assertTrue(dec(time_arr, self.delta)) - self.assertAlmostEqual( - np.abs(np.max(steady - transient[-1, :])), 0.0, places=4 - ) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=4) def test_ext_3d(self): steady = af.ext_thiem_3d( @@ -127,9 +123,7 @@ def test_ext_3d(self): self.assertTrue(inc(rad_arr, self.delta)) for time_arr in transient.T: self.assertTrue(dec(time_arr, self.delta)) - self.assertAlmostEqual( - np.abs(np.max(steady - transient[-1, :])), 0.0, places=4 - ) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=4) def test_neuman(self): steady = af.neuman2004_steady( @@ -157,9 +151,7 @@ def test_neuman(self): self.assertTrue(inc(rad_arr, self.delta)) for time_arr in transient.T: self.assertTrue(dec(time_arr, self.delta)) - self.assertAlmostEqual( - np.abs(np.max(steady - transient[-1, :])), 0.0, places=4 - ) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=4) def test_tpl(self): steady = af.ext_thiem_tpl( @@ -193,9 +185,7 @@ def test_tpl(self): self.assertTrue(inc(rad_arr, self.delta)) for time_arr in transient.T: self.assertTrue(dec(time_arr, self.delta)) - self.assertAlmostEqual( - np.abs(np.max(steady - transient[-1, :])), 0.0, places=3 - ) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=3) def test_tpl_3d(self): steady = af.ext_thiem_tpl_3d( @@ -229,9 +219,75 @@ def test_tpl_3d(self): self.assertTrue(inc(rad_arr, self.delta)) for time_arr in transient.T: self.assertTrue(dec(time_arr, self.delta)) - self.assertAlmostEqual( - np.abs(np.max(steady - transient[-1, :])), 0.0, places=2 + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=2) + + def test_int(self): + steady = af.ext_thiem_int( + rad=self.rad, + r_ref=self.r_ref, + cond_gmean=self.t_gmean, + var=self.var, + len_scale=self.len_scale, + lat_ext=self.lat_ext, + roughness=self.hurst, + dim=self.frac_dim, + rate=self.rate, + h_ref=self.h_ref, + ) + transient = af.ext_theis_int( + time=self.time, + rad=self.rad, + storage=self.storage, + cond_gmean=self.t_gmean, + var=self.var, + len_scale=self.len_scale, + lat_ext=self.lat_ext, + roughness=self.hurst, + dim=self.frac_dim, + rate=self.rate, + r_bound=self.r_ref, + h_bound=self.h_ref, + ) + self.assertTrue(inc(steady, self.delta)) + for rad_arr in transient: + self.assertTrue(inc(rad_arr, self.delta)) + for time_arr in transient.T: + self.assertTrue(dec(time_arr, self.delta)) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=3) + + def test_int_3d(self): + steady = af.ext_thiem_int_3d( + rad=self.rad, + r_ref=self.r_ref, + cond_gmean=self.t_gmean, + var=self.var, + len_scale=self.len_scale, + lat_ext=self.lat_ext, + roughness=self.hurst, + anis=self.anis, + rate=self.rate, + h_ref=self.h_ref, + ) + transient = af.ext_theis_int_3d( + time=self.time, + rad=self.rad, + storage=self.storage, + cond_gmean=self.t_gmean, + var=self.var, + len_scale=self.len_scale, + lat_ext=self.lat_ext, + roughness=self.hurst, + anis=self.anis, + rate=self.rate, + r_bound=self.r_ref, + h_bound=self.h_ref, ) + self.assertTrue(inc(steady, self.delta)) + for rad_arr in transient: + self.assertTrue(inc(rad_arr, self.delta)) + for time_arr in transient.T: + self.assertTrue(dec(time_arr, self.delta)) + self.assertAlmostEqual(np.abs(np.max(steady - transient[-1, :])), 0.0, places=2) # plt.plot(self.rad, steady) # for rad_arr in transient: