diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 2a2e539..d4033ef 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -1,2 +1,6 @@ # Migrate code style to black -998b86fc128e35f1506e3d4facb169c558679bde \ No newline at end of file +998b86fc128e35f1506e3d4facb169c558679bde + +# Migrate code style to ruff +a11f92ff2732d7fdc548559f15f009e4a3d5ea7a +571e4b5af45bd1ca6f025dd883b9ea800acb4aef \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b2c570f..7099b4d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -30,12 +30,10 @@ jobs: python -m pip install --upgrade pip pip install -e .[ci] - - name: Lint with flake8 + - name: Lint with ruff run: | - # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=80 --statistics + ruff check . + ruff format --check - name: Run tests run: | diff --git a/examples/timml_notebook0_sol.py b/examples/timml_notebook0_sol.py index d733184..2b6c999 100644 --- a/examples/timml_notebook0_sol.py +++ b/examples/timml_notebook0_sol.py @@ -1,34 +1,32 @@ - # coding: utf-8 # # TimML Notebook 0 # ## Single layer flow -# Consider uniform flow from East to West. The gradient is 0.001. The hydraulic conductivity is $k=10$ m/d. The aquifer bottom and top are at 0 m and 10 m. The head at $x=-1000$ m and $y=0$ is fixed at 41 m. +# Consider uniform flow from East to West. The gradient is 0.001. The hydraulic conductivity is $k=10$ m/d. The aquifer bottom and top are at 0 m and 10 m. The head at $x=-1000$ m and $y=0$ is fixed at 41 m. # In[1]: +import numpy as np -from timml import * -from pylab import * - +import timml # In[2]: -ml = ModelMaq(kaq=10, z=[10, 0]) +ml = timml.ModelMaq(kaq=10, z=[10, 0]) # In[3]: -rf = Constant(ml, xr=-1000, yr=0, hr=41) +rf = timml.Constant(ml, xr=-1000, yr=0, hr=41) # In[4]: -uf = Uflow(ml, slope=0.001, angle=0) +uf = timml.Uflow(ml, slope=0.001, angle=0) # In[5]: @@ -40,18 +38,37 @@ # In[6]: -ml.contour(x1=-1200, x2=200, nx=50, y1=-500, y2=500, ny=50, - levels=10, labels=True, decimals=2, legend=True) +ml.contour( + x1=-1200, + x2=200, + nx=50, + y1=-500, + y2=500, + ny=50, + levels=10, + labels=True, + decimals=2, + legend=True, +) -# The default contour levels are not what we want for this example, so let's specify the levels +# The default contour levels are not what we want for this example, so let's specify the levels # to go from 39 to 42 with steps of 0.1 (not all those levels are present in the current window). # In[7]: -ml.contour(x1=-1200, x2=200, nx=50, y1=-500, y2=500, ny=50, - levels=arange(39, 42, 0.1), labels=True, decimals=1) +ml.contour( + x1=-1200, + x2=200, + nx=50, + y1=-500, + y2=500, + ny=50, + levels=np.arange(39, 42, 0.1), + labels=True, + decimals=1, +) # A well is located at $(x,y)=(-400,0)$ with a discharge $Q=50$ m$^3$ and a radius of 0.2 m. @@ -59,17 +76,17 @@ # In[8]: -w = Well(ml, xw=-400, yw=0, Qw=50., rw=0.2) +w = timml.Well(ml, xw=-400, yw=0, Qw=50.0, rw=0.2) -# After the well is added (or any other elements), the model needs to be solved again. A contour plot is created and a 10 strace line are added. The stepsize is given in meters and represents the largest space step that is taken, but it is reduced when certain accuracy constraints are not met. Note that, after running the code cell below, for each trace line it is printed to the screen what the reason was that the traceline was aborted. In this case it was either because the trace line reached a well or reached the maximum number of steps (the default is 100 steps, but this can be changed by specifying the `nstepmax` keyword). +# After the well is added (or any other elements), the model needs to be solved again. A contour plot is created and a 10 strace line are added. The stepsize is given in meters and represents the largest space step that is taken, but it is reduced when certain accuracy constraints are not met. Note that, after running the code cell below, for each trace line it is printed to the screen what the reason was that the traceline was aborted. In this case it was either because the trace line reached a well or reached the maximum number of steps (the default is 100 steps, but this can be changed by specifying the `nstepmax` keyword). # In[9]: ml.solve() ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1)) -ml.tracelines(-800 * ones(1), -200 * ones(1), zeros(1), hstepmax=20) +ml.tracelines(-800 * np.ones(1), -200 * np.ones(1), np.zeros(1), hstepmax=20) # ### Exercise a @@ -79,7 +96,7 @@ ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1)) -ml.tracelines(-800 * ones(10), linspace(-500, 500, 10), zeros(10), hstepmax=20) +ml.tracelines(-800 * np.ones(10), np.linspace(-500, 500, 10), np.zeros(10), hstepmax=20) # ### Exercise b @@ -88,14 +105,14 @@ # In[11]: -ml = ModelMaq(kaq=10, z=[10, 0]) -rf = Constant(ml, xr=-1000, yr=0, hr=41) -uf = Uflow(ml, slope=0.001, angle=0) -w = Well(ml, xw=-400, yw=0, Qw=200, rw=0.2) +ml = timml.ModelMaq(kaq=10, z=[10, 0]) +rf = timml.Constant(ml, xr=-1000, yr=0, hr=41) +uf = timml.Uflow(ml, slope=0.001, angle=0) +w = timml.Well(ml, xw=-400, yw=0, Qw=200, rw=0.2) ml.solve() ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1)) -ml.tracelines(-800 * ones(10), linspace(-500, 500, 10), zeros(10), hstepmax=20) -print(('head at well:', w.headinside())) +ml.tracelines(-800 * np.ones(10), np.linspace(-500, 500, 10), np.zeros(10), hstepmax=20) +print(("head at well:", w.headinside())) # ### Add a river @@ -104,33 +121,37 @@ # In[12]: -ml = ModelMaq(kaq=10, z=[10, 0]) -rf = Constant(ml, xr=-1000, yr=0, hr=41) -uf = Uflow(ml, slope=0.001, angle=0) -w = Well(ml, xw=-400, yw=0, Qw=200, rw=0.2) -ls1 = HeadLineSink(ml, 0, -500, 0, 500, 40) +ml = timml.ModelMaq(kaq=10, z=[10, 0]) +rf = timml.Constant(ml, xr=-1000, yr=0, hr=41) +uf = timml.Uflow(ml, slope=0.001, angle=0) +w = timml.Well(ml, xw=-400, yw=0, Qw=200, rw=0.2) +ls1 = timml.HeadLineSink(ml, 0, -500, 0, 500, 40) ml.solve() -ml.contour(-1000, 100, 50, -500, 500, 50, levels=arange(39, 42, 0.1)) -print(('head at well:', w.headinside())) +ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1)) +print(("head at well:", w.headinside())) # ### Exercise c -# Simulate the river with 20 line-sinks from $y=-800$ to $y=800$. +# Simulate the river with 20 line-sinks from $y=-800$ to $y=800$. # In[13]: -ml = ModelMaq(kaq=10, z=[10, 0]) -rf = Constant(ml, xr=-1000, yr=0, hr=41) -uf = Uflow(ml, slope=0.001, angle=0) -w = Well(ml, xw=-400, yw=0, Qw=200, rw=0.2) -xls = zeros(21) -yls = linspace(-800, 800, 21) -ls = HeadLineSinkString(ml, xy=list(zip(xls, yls)), hls=40, layers=0) +ml = timml.ModelMaq(kaq=10, z=[10, 0]) +rf = timml.Constant(ml, xr=-1000, yr=0, hr=41) +uf = timml.Uflow(ml, slope=0.001, angle=0) +w = timml.Well(ml, xw=-400, yw=0, Qw=200, rw=0.2) +xls = np.zeros(21) +yls = np.linspace(-800, 800, 21) +ls = timml.HeadLineSinkString(ml, xy=list(zip(xls, yls)), hls=40, layers=0) ml.solve() -ml.contour(-1000, 100, 50, -500, 500, 50, levels=arange(39, 42, 0.1)) -ml.tracelines(-800 * ones(10), linspace(-500, 500, 10), zeros(10), hstepmax=20, color='b') -ml.tracelines(-0.01 * ones(5), linspace(-150, 150, 5), zeros(5), hstepmax=20, color='r') +ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1)) +ml.tracelines( + -800 * np.ones(10), np.linspace(-500, 500, 10), np.zeros(10), hstepmax=20, color="b" +) +ml.tracelines( + -0.01 * np.ones(5), np.linspace(-150, 150, 5), np.zeros(5), hstepmax=20, color="r" +) # ### Capture zone @@ -139,8 +160,8 @@ # In[14]: -ml.contour(-1000, 100, 50, -500, 500, 50, levels=arange(39, 42, 0.1), layers=0) -w.plotcapzone(hstepmax=20, nt=20, zstart=0, tmax=5 * 365.25, color='b') +ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1), layers=0) +w.plotcapzone(hstepmax=20, nt=20, zstart=0, tmax=5 * 365.25, color="b") # ### Exercise d @@ -149,12 +170,8 @@ # In[15]: -ml.contour(-1000, 100, 50, -500, 500, 50, levels=arange(39, 42, 0.1), layers=0) -w.plotcapzone(hstepmax=20, nt=20, zstart=0, tmax=20 * 365.25, color='b') +ml.contour(-1000, 100, 50, -500, 500, 50, levels=np.arange(39, 42, 0.1), layers=0) +w.plotcapzone(hstepmax=20, nt=20, zstart=0, tmax=20 * 365.25, color="b") # In[ ]: - - - - diff --git a/examples/workshop_nhv/vlaketunnel_functions.py b/examples/workshop_nhv/vlaketunnel_functions.py index a640ab0..eeff2c0 100644 --- a/examples/workshop_nhv/vlaketunnel_functions.py +++ b/examples/workshop_nhv/vlaketunnel_functions.py @@ -1,13 +1,14 @@ # -*- coding: utf-8 -*- -import numpy as np import matplotlib.pyplot as plt +import numpy as np + import timml as tml -def create_model(kaq=[0.1, 5.0, 15.0, 5.0], c=[1000.0, 2.0, 2.0, 2.0],hstar=0, c_channel_bot=30, - do_plot=True, df_dh=None): + +def create_model(kaq=None, c=None, hstar=0, c_channel_bot=30, do_plot=True, df_dh=None): """ - Create a TimML model for Vlaketunnel case + Create a TimML model for Vlaketunnel case. Parameters ---------- @@ -30,8 +31,11 @@ def create_model(kaq=[0.1, 5.0, 15.0, 5.0], c=[1000.0, 2.0, 2.0, 2.0],hstar=0, c The model """ - # create model + if c is None: + c = [1000.0, 2.0, 2.0, 2.0] + if kaq is None: + kaq = [0.1, 5.0, 15.0, 5.0] ml = tml.ModelMaq( kaq=kaq, z=[1.0, -3.0, -7.0, -7.0, -14.0, -14.0, -30.0, -30.0, -40.0], @@ -42,20 +46,30 @@ def create_model(kaq=[0.1, 5.0, 15.0, 5.0], c=[1000.0, 2.0, 2.0, 2.0],hstar=0, c ) # add dewatering - dewatering_east_xys = [[59224, 387382], [59359, 387375], [59360, 387311], [59234, 387298], ] - q_east_total = 325*24 - - q_west_total = 75*24 - dewatering_west_xys = [[58781, 387375], [58785, 387307],] - - for dewatering_xys, q_total in zip([dewatering_east_xys, dewatering_west_xys], [q_east_total, q_west_total]): + dewatering_east_xys = [ + [59224, 387382], + [59359, 387375], + [59360, 387311], + [59234, 387298], + ] + q_east_total = 325 * 24 + + q_west_total = 75 * 24 + dewatering_west_xys = [ + [58781, 387375], + [58785, 387307], + ] + + for dewatering_xys, q_total in zip( + [dewatering_east_xys, dewatering_west_xys], [q_east_total, q_west_total] + ): # loop over both dewatering locations for dewatering_xy in dewatering_xys: # loop over the modelled wells, in pratice a lot of more wells are used. Current model has focus on regional effect, therefore limited number of wells are considered sufficient - dewatering_east = tml.Well( + tml.Well( xw=dewatering_xy[0], yw=dewatering_xy[1], - Qw=q_total/len(dewatering_xys), + Qw=q_total / len(dewatering_xys), rw=0.5, res=1.0, layers=1, @@ -66,7 +80,7 @@ def create_model(kaq=[0.1, 5.0, 15.0, 5.0], c=[1000.0, 2.0, 2.0, 2.0],hstar=0, c c_channel = ml.aq.c.copy() c_channel[0] = c_channel_bot - channel_0 = tml.PolygonInhomMaq( + tml.PolygonInhomMaq( kaq=ml.aq.kaq, z=ml.aq.z, c=c_channel, @@ -74,21 +88,32 @@ def create_model(kaq=[0.1, 5.0, 15.0, 5.0], c=[1000.0, 2.0, 2.0, 2.0],hstar=0, c npor=[None, None, None, None, None, None, None, None], hstar=0.0, # compared to QGIS-Tim export the channel is extended to the north in order to cover the northern observation wells better - xy= [ [58921, 390500], [59065, 390500], [59110, 387996], [59146, 387447], [59263, 386809], [59317, 386260], [59110, 386251], [58966, 386863], [58921, 388617], ], + xy=[ + [58921, 390500], + [59065, 390500], + [59110, 387996], + [59146, 387447], + [59263, 386809], + [59317, 386260], + [59110, 386251], + [58966, 386863], + [58921, 388617], + ], order=4, ndeg=6, model=ml, ) ml.solve() - + if do_plot and (df_dh is not None): plot_model_results(ml, df_dh) - + return ml + def plot_model_input(ml): """ - Plot model input in schematic section + Plot model input in schematic section. Parameters ---------- @@ -101,43 +126,55 @@ def plot_model_input(ml): """ # some plotting constants - xmin=-1 - xchannel=-0.25 - xhinter=-0.2 - xmax=1 - zaqmid = np.mean([ml.aq.zaqtop,ml.aq.zaqbot],axis=0) + xmin = -1 + xchannel = -0.25 + xhinter = -0.2 + xmax = 1 + zaqmid = np.mean([ml.aq.zaqtop, ml.aq.zaqbot], axis=0) # plot layers - plt.hlines(y=ml.aq.zlltop,xmin=xmin,xmax=xmax,color='darkgray') - plt.hlines(y=ml.aq.zaqbot,xmin=xmin,xmax=xmax,color='darkgray') + plt.hlines(y=ml.aq.zlltop, xmin=xmin, xmax=xmax, color="darkgray") + plt.hlines(y=ml.aq.zaqbot, xmin=xmin, xmax=xmax, color="darkgray") # plot kh for kh, z in zip(ml.aq.kaq, zaqmid): - plt.annotate(f'kh={kh:0.1f}m/d',(0,z),ha='center') + plt.annotate(f"kh={kh:0.1f}m/d", (0, z), ha="center") # plot c for c, z in zip(ml.aq.c, ml.aq.zaqtop): - plt.annotate(f'c={c:0.1f}d',(0.5,z),ha='center',va='center') + plt.annotate(f"c={c:0.1f}d", (0.5, z), ha="center", va="center") # plot channel - plt.plot([xmin,xchannel],[ml.aq.inhomlist[0].hstar]*2,color='blue') - plt.annotate(f'h_ch={ml.aq.inhomlist[0].hstar:0.1f}',(xchannel,ml.aq.inhomlist[0].hstar),ha='right',va='bottom') - plt.annotate(f'c_ch={ml.aq.inhomlist[0].c[0]:0.1f}',(xchannel,ml.aq.zaqtop[0]),ha='right',va='bottom') + plt.plot([xmin, xchannel], [ml.aq.inhomlist[0].hstar] * 2, color="blue") + plt.annotate( + f"h_ch={ml.aq.inhomlist[0].hstar:0.1f}", + (xchannel, ml.aq.inhomlist[0].hstar), + ha="right", + va="bottom", + ) + plt.annotate( + f"c_ch={ml.aq.inhomlist[0].c[0]:0.1f}", + (xchannel, ml.aq.zaqtop[0]), + ha="right", + va="bottom", + ) # plot hinterland - plt.plot([xhinter,xmax],[ml.aq.hstar]*2,color='darkblue') - plt.annotate(f'h_polder={ml.aq.hstar:0.1f}',(xhinter,ml.aq.hstar),ha='left',va='bottom') + plt.plot([xhinter, xmax], [ml.aq.hstar] * 2, color="darkblue") + plt.annotate( + f"h_polder={ml.aq.hstar:0.1f}", (xhinter, ml.aq.hstar), ha="left", va="bottom" + ) plt.xlim([xmin, xmax]) def plot_model_results(ml, df_dh): """ - Plot results of TimML model of Vlaketunnel case + Plot results of TimML model of Vlaketunnel case. Parameters ---------- - ml : timml Model, + ml : timml Model, The model. - df_dh : pd.DataFrame + df_dh : pd.DataFrame Observed drawdowns Returns @@ -145,37 +182,50 @@ def plot_model_results(ml, df_dh): None. """ - # contour plot plt.subplot(221) - ml.contour(win=[57000, 60000, 386900, 389100], ngr=50, layers=1, - levels=[-5,-2,-1,-0.5,-0.1], labels=True, decimals=2, legend=False, newfig=False); + ml.contour( + win=[57000, 60000, 386900, 389100], + ngr=50, + layers=1, + levels=[-5, -2, -1, -0.5, -0.1], + labels=True, + decimals=2, + legend=False, + newfig=False, + ) plt.scatter(df_dh.x, df_dh.y, 20, c=df_dh.color) - for index, row in df_dh.iterrows(): - plt.annotate(f'{row.dh_obs:0.2f}', (row.x, row.y),ha=row.ha,va=row.va) - plt.title('contours in layer 1'); - + for _, row in df_dh.iterrows(): + plt.annotate(f"{row.dh_obs:0.2f}", (row.x, row.y), ha=row.ha, va=row.va) + plt.title("contours in layer 1") + # plot model input plt.subplot(222) plot_model_input(ml) - + for plotid in (223, 224): plt.subplot(plotid) if plotid == 223: # first plot, get model results - df_dh['ml_layer'] = None - df_dh['dh_calc'] = None + df_dh["ml_layer"] = None + df_dh["dh_calc"] = None for index, row in df_dh.iterrows(): - df_dh.loc[index,'ml_layer'] = np.where(ml.aq.zaqtop > row.screen_top)[0][-1] - df_dh.loc[index,'dh_calc'] = ml.headalongline(row.x, row.y, row.ml_layer)[0][0] + df_dh.loc[index, "ml_layer"] = np.where(ml.aq.zaqtop > row.screen_top)[ + 0 + ][-1] + df_dh.loc[index, "dh_calc"] = ml.headalongline( + row.x, row.y, row.ml_layer + )[0][0] # plot all model results plot_df = df_dh else: # second plot, only plot outside dewatering area plot_df = df_dh.loc[df_dh.r > 100] - plt.scatter(plot_df.r, plot_df.dh_obs, 50, c=plot_df.color, alpha=0.3, label='observed') - plt.scatter(plot_df.r, plot_df.dh_calc, 40, marker='+', label='modelled') + plt.scatter( + plot_df.r, plot_df.dh_obs, 50, c=plot_df.color, alpha=0.3, label="observed" + ) + plt.scatter(plot_df.r, plot_df.dh_calc, 40, marker="+", label="modelled") plt.legend() - plt.title('heads from screened modellayer'); - plt.grid() \ No newline at end of file + plt.title("heads from screened modellayer") + plt.grid() diff --git a/pyproject.toml b/pyproject.toml index 4b882e2..bb36c2a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,7 +54,7 @@ ci = [ "pytest>=4.6", "pytest-cov", "jupyter>=1.0.0", - "flake8", + "ruff", "coveralls", "shapely", ] @@ -74,10 +74,30 @@ packages = ["timml", "timml.besselaesnumba"] [tool.setuptools.dynamic] version = { attr = "timml.version.__version__" } -[tool.black] -line-length = 88 +[tool.ruff.lint] +# See: https://docs.astral.sh/ruff/rules/ +select = [ + "C4", # flake8-comprehensions + "E", # pycodestyle + "F", # pyflakes + "I", # isort + "PT", # pytest-style + "D", # pydocstyle + "B", # flake8-bugbear + "NPY", # numpy +] +ignore = [ + "D100", # Missing docstring in public module + "E501", # Line too long + "D101", # Missing docstring in public class + "D102", # Missing docstring in public method + "D103", # Missing docstring in public function + "D104", # Missing docstring in public package + "D105", # Missing docstring in magic method + "D205", # 1 blank line required between summary line and description + "D401", # First line of docstring should be in imperative mood + "E741", # Ambiguous variable name (such as "l") +] -[tool.isort] -profile = "black" -src_paths = ["timml"] -line_length = 88 +[tool.ruff.lint.pydocstyle] +convention = "numpy" \ No newline at end of file diff --git a/tests/test_besselaes.py b/tests/test_besselaes.py index a9edca0..0093512 100644 --- a/tests/test_besselaes.py +++ b/tests/test_besselaes.py @@ -1,5 +1,4 @@ import numpy as np -import pytest from numpy.testing import assert_allclose from timml.besselaesnumba import besselaesnumba as besselaesnew diff --git a/tests/test_import.py b/tests/test_import.py index ff8b028..dbc7b77 100644 --- a/tests/test_import.py +++ b/tests/test_import.py @@ -1,6 +1,8 @@ def test_import(): import timml + print(timml.__version__) + if __name__ == "__main__": test_import() diff --git a/tests/test_notebooks.py b/tests/test_notebooks.py index ba5d1a3..ec7108b 100644 --- a/tests/test_notebooks.py +++ b/tests/test_notebooks.py @@ -28,24 +28,21 @@ def get_notebooks(): def get_jupyter_kernel(): - try: - jklcmd = ("jupyter", "kernelspec", "list") - b = subprocess.Popen( - jklcmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT - ).communicate()[0] - if isinstance(b, bytes): - b = b.decode("utf-8") - print(b) - for line in b.splitlines(): - if "python" in line: - kernel = line.split()[0] - except: - kernel = None + jklcmd = ("jupyter", "kernelspec", "list") + b = subprocess.Popen( + jklcmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT + ).communicate()[0] + if isinstance(b, bytes): + b = b.decode("utf-8") + print(b) + for line in b.splitlines(): + if "python" in line: + kernel = line.split()[0] return kernel -@pytest.mark.notebooks +@pytest.mark.notebooks() @pytest.mark.parametrize("pth", get_notebooks()) def test_notebook(pth): kernel = get_jupyter_kernel() diff --git a/timml/__init__.py b/timml/__init__.py index 054f3b0..b60dbc9 100644 --- a/timml/__init__.py +++ b/timml/__init__.py @@ -1,4 +1,5 @@ -"""Copyright (C), 2015, Mark Bakker. +"""Copyright (C), 2015, Mark Bakker. + Mark Bakker, Delft University of Technology mark dot bakker at tudelft dot nl @@ -10,12 +11,12 @@ # --version number __name__ = "timml" __author__ = "Mark Bakker" -from . import bessel +from timml import bessel # Import all classes and functions -from .circareasink import CircAreaSink -from .constant import Constant, ConstantStar -from .inhomogeneity import ( +from timml.circareasink import CircAreaSink +from timml.constant import Constant, ConstantStar +from timml.inhomogeneity import ( BuildingPit3D, BuildingPitMaq, LeakyBuildingPit3D, @@ -23,15 +24,15 @@ PolygonInhom3D, PolygonInhomMaq, ) -from .inhomogeneity1d import StripInhom3D, StripInhomMaq -from .linedoublet import ( +from timml.inhomogeneity1d import StripInhom3D, StripInhomMaq +from timml.linedoublet import ( ImpLineDoublet, ImpLineDoubletString, LeakyLineDoublet, LeakyLineDoubletString, ) -from .linedoublet1d import ImpLineDoublet1D, LeakyLineDoublet1D -from .linesink import ( +from timml.linedoublet1d import ImpLineDoublet1D, LeakyLineDoublet1D +from timml.linesink import ( HeadLineSink, HeadLineSinkContainer, HeadLineSinkString, @@ -40,15 +41,54 @@ LineSinkDitch, LineSinkDitchString, ) -from .linesink1d import HeadLineSink1D, LineSink1D -from .model import Model, Model3D, ModelMaq -from .stripareasink import StripAreaSink -from .trace import timtraceline, timtracelines -from .uflow import Uflow -from .version import __version__ -from .well import HeadWell, LargeDiameterWell, Well, WellBase +from timml.linesink1d import HeadLineSink1D, LineSink1D +from timml.model import Model, Model3D, ModelMaq +from timml.stripareasink import StripAreaSink +from timml.trace import timtraceline, timtracelines +from timml.uflow import Uflow +from timml.version import __version__ +from timml.well import HeadWell, LargeDiameterWell, Well, WellBase -__all__ = [s for s in dir() if not s.startswith("_")] +__all__ = [ + "CircAreaSink", + "Constant", + "ConstantStar", + "BuildingPit3D", + "BuildingPitMaq", + "LeakyBuildingPit3D", + "LeakyBuildingPitMaq", + "PolygonInhom3D", + "PolygonInhomMaq", + "StripInhom3D", + "StripInhomMaq", + "ImpLineDoublet", + "ImpLineDoubletString", + "LeakyLineDoublet", + "LeakyLineDoubletString", + "ImpLineDoublet1D", + "LeakyLineDoublet1D", + "HeadLineSink", + "HeadLineSinkContainer", + "HeadLineSinkString", + "HeadLineSinkZero", + "LineSinkBase", + "LineSinkDitch", + "LineSinkDitchString", + "HeadLineSink1D", + "LineSink1D", + "Model", + "Model3D", + "ModelMaq", + "StripAreaSink", + "timtraceline", + "timtracelines", + "Uflow", + "__version__", + "HeadWell", + "LargeDiameterWell", + "Well", + "WellBase", +] # default bessel module is numba bessel.set_bessel_method(method="numba") diff --git a/timml/aquifer.py b/timml/aquifer.py index 05dabc5..3ca9c43 100644 --- a/timml/aquifer.py +++ b/timml/aquifer.py @@ -2,7 +2,6 @@ import numpy as np -from .aquifer_parameters import param_maq from .constant import ConstantStar @@ -33,9 +32,9 @@ def __init__(self, model, kaq, c, z, npor, ltype): self.layernumber[self.ltype == "a"] = np.arange(self.naq) self.layernumber[self.ltype == "l"] = np.arange(self.nlayers - self.naq) if self.ltype[0] == "a": - self.layernumber[ - self.ltype == "l" - ] += 1 # first leaky layer below first aquifer layer + self.layernumber[self.ltype == "l"] += ( + 1 # first leaky layer below first aquifer layer + ) self.zaqtop = self.z[:-1][self.ltype == "a"] self.zaqbot = self.z[1:][self.ltype == "a"] self.Haq = self.zaqtop - self.zaqbot diff --git a/timml/bessel.py b/timml/bessel.py index 3473e50..edee416 100644 --- a/timml/bessel.py +++ b/timml/bessel.py @@ -10,7 +10,11 @@ def set_bessel_method(method="numba"): bessel = besselaesnew.besselaesnew bessel.initialize() except ImportError: - warn("Cannot import compiled fortran bessel module! Defaulting to numba!") + warn( + "Cannot import compiled fortran bessel module! Defaulting to numba!", + category=ImportWarning, + stacklevel=1, + ) bessel = import_module("timml.besselaesnumba.besselaesnumba") elif method == "numba": bessel = import_module("timml.besselaesnumba.besselaesnumba") diff --git a/timml/besselaesnumba/besselaesnumba.py b/timml/besselaesnumba/besselaesnumba.py index 78997b1..8993084 100644 --- a/timml/besselaesnumba/besselaesnumba.py +++ b/timml/besselaesnumba/besselaesnumba.py @@ -110,11 +110,11 @@ def potbeslsho(x, y, z1, z2, labda, order, ilap, naq): ilap: equals 1 when first value is Laplace line-sink and first labda equals zero naq: Number of aquifers rv(naq): Array to store return value (must be pre-allocated) + Returns - -------- + ------- rv(naq): Potentials. First spot is Laplace value if ilap=1 """ - rv = np.zeros(naq) # lstype = 1 means line-sink diff --git a/timml/besselaesnumba/besselaesnumba_old.py b/timml/besselaesnumba/besselaesnumba_old.py index d458c87..2343523 100644 --- a/timml/besselaesnumba/besselaesnumba_old.py +++ b/timml/besselaesnumba/besselaesnumba_old.py @@ -99,11 +99,11 @@ def potbeslsho(x, y, z1, z2, labda, order, ilap, naq): ilap: equals 1 when first value is Laplace line-sink and first labda equals zero naq: Number of aquifers rv(naq): Array to store return value (must be pre-allocated) + Returns - -------- + ------- rv(naq): Potentials. Fist spot is Laplace value if ilap=1 """ - rv = np.zeros(naq) # lstype = 1 means line-sink diff --git a/timml/circareasink.py b/timml/circareasink.py index 3cee255..2786d4d 100644 --- a/timml/circareasink.py +++ b/timml/circareasink.py @@ -246,8 +246,5 @@ def changetrace( u = u1 * (1.0 + eps) # Go just beyond circle else: u = u2 * (1.0 + eps) # Go just beyond circle - xn = x1 + u * (x2 - x1) - yn = y1 + u * (y2 - y1) - zn = xyzt1[2] + u * (xyzt2[2] - xyzt1[2]) xyztnew = xyzt1 + u * (xyzt2 - xyzt1) return changed, terminate, xyztnew, message diff --git a/timml/circinhom.py b/timml/circinhom.py index 813d5dc..595a473 100644 --- a/timml/circinhom.py +++ b/timml/circinhom.py @@ -4,9 +4,9 @@ TimML.py file for more details. (c) Mark Bakker, 2002-2007 """ +import numpy as np import scipy.special from element import Element -import numpy as np class CircleInhom(Element): @@ -133,7 +133,7 @@ def potentialInfluence(self, aq, x, y): def potentialInfluenceInLayer(self, aq, pylayer, x, y): """Returns PotentialInfluence function in aquifer aq in pylayer as array (1 - value per parameter) + value per parameter). Notes ----- @@ -152,7 +152,8 @@ def potentialInfluenceInLayer(self, aq, pylayer, x, y): def potentialInfluenceAllLayers(self, aq, pylayer, x, y): """Returns PotentialInfluence function in aquifer aq in all layers as an - array.""" + array. + """ potInf = self.potentialInfluence(aq, x, y) rv = np.zeros((aq.Naquifers, 0), "d") for p in potInf: @@ -165,7 +166,8 @@ def potentialInfluenceAllLayers(self, aq, pylayer, x, y): def potentialInfluenceSpecLayers(self, aq, pylayer, x, y): """Returns PotentialInfluence function in aquifer aq in all layers as an - array.""" + array. + """ potInf = self.potentialInfluence(aq, x, y) pylen = len(pylayer) rv = np.zeros((pylen, 0), "d") diff --git a/timml/element.py b/timml/element.py index 62c6726..c4156d0 100644 --- a/timml/element.py +++ b/timml/element.py @@ -27,7 +27,7 @@ def initialize(self): pass def potinf(self, x, y, aq=None): - """Returns array of size (nparam, naq)""" + """Returns array of size (nparam, naq).""" raise Exception("Must overload Element.potinf()") def potential(self, x, y, aq=None): @@ -37,7 +37,8 @@ def potential(self, x, y, aq=None): def potinflayers(self, x, y, layers, aq=None): """Returns array of size (len(layers),nparam) only used in building - equations.""" + equations. + """ if aq is None: aq = self.model.aq.find_aquifer_data(x, y) pot = self.potinf(x, y, aq) # nparam rows, naq cols @@ -54,18 +55,19 @@ def potentiallayers(self, x, y, layers, aq=None): return pot[layers] def disvecinf(self, x, y, aq=None): - """Returns array of size (2, nparam, naq)""" + """Returns array of size (2, nparam, naq).""" raise Exception("Must overload Element.disvecinf()") def disvec(self, x, y, aq=None): - """Returns array of size (2, nparam, naq)""" + """Returns array of size (2, nparam, naq).""" if aq is None: aq = self.model.aq.find_aquifer_data(x, y) return np.sum(self.parameters * self.disvecinf(x, y, aq), 1) def disvecinflayers(self, x, y, layers, aq=None): """Returns two arrays of size (len(layers),nparam) only used in building - equations.""" + equations. + """ if aq is None: aq = self.model.aq.find_aquifer_data(x, y) qxqy = self.disvecinf(x, y, aq) # nparam rows, naq cols diff --git a/timml/equation.py b/timml/equation.py index ba9f557..374f1e7 100644 --- a/timml/equation.py +++ b/timml/equation.py @@ -21,9 +21,9 @@ def equation(self): ieq = 0 for e in self.model.elementlist: if e.nunknowns > 0: - mat[ - istart : istart + self.nlayers, ieq : ieq + e.nunknowns - ] = e.potinflayers(self.xc[icp], self.yc[icp], self.layers) + mat[istart : istart + self.nlayers, ieq : ieq + e.nunknowns] = ( + e.potinflayers(self.xc[icp], self.yc[icp], self.layers) + ) if e == self: mat[ istart : istart + self.nlayers, ieq : ieq + e.nunknowns @@ -100,9 +100,9 @@ def equation(self): ieq = 0 for e in self.model.elementlist: if e.nunknowns > 0: - mat[ - istart : istart + self.nlayers, ieq : ieq + e.nunknowns - ] = e.potinflayers(self.xc[icp], self.yc[icp], self.layers) + mat[istart : istart + self.nlayers, ieq : ieq + e.nunknowns] = ( + e.potinflayers(self.xc[icp], self.yc[icp], self.layers) + ) ieq += e.nunknowns else: rhs[istart : istart + self.nlayers] -= e.potentiallayers( diff --git a/timml/inhomogeneity.py b/timml/inhomogeneity.py index 8e4afd3..a27b929 100644 --- a/timml/inhomogeneity.py +++ b/timml/inhomogeneity.py @@ -79,7 +79,7 @@ def create_elements(self): self.zcout[i].real, self.zcout[i].imag ) if (aqout == self.model.aq) or (aqout.inhom_number > self.inhom_number): - ls = IntHeadDiffLineSink( + IntHeadDiffLineSink( self.model, x1=self.x[i], y1=self.y[i], @@ -93,7 +93,7 @@ def create_elements(self): aqin=aqin, aqout=aqout, ) - ls = IntFluxDiffLineSink( + IntFluxDiffLineSink( self.model, x1=self.x[i], y1=self.y[i], @@ -171,8 +171,8 @@ def __init__( model, xy, kaq=1, - z=[1, 0], - c=[], + z=None, + c=None, npor=0.3, topboundary="conf", hstar=None, @@ -180,6 +180,10 @@ def __init__( order=3, ndeg=3, ): + if c is None: + c = [] + if z is None: + z = [1, 0] if N is not None: assert ( topboundary[:4] == "conf" @@ -251,7 +255,7 @@ def __init__( model, xy, kaq=1, - z=[1, 0], + z=None, kzoverkh=1, npor=0.3, topboundary="conf", @@ -262,6 +266,8 @@ def __init__( order=3, ndeg=3, ): + if z is None: + z = [1, 0] if N is not None: assert ( topboundary[:4] == "conf" @@ -311,7 +317,7 @@ def __init__( npor=0.3, order=3, ndeg=3, - layers=[0], + layers=None, ): """Element to simulate a building pit with an impermeable wall. @@ -355,6 +361,8 @@ def __init__( layers: list or np.array layers in which impermeable wall is present. """ + if layers is None: + layers = [0] AquiferData.__init__(self, model, kaq, c, z, npor, ltype) self.order = order self.ndeg = ndeg @@ -501,14 +509,14 @@ def __init__( model, xy, kaq=1.0, - c=[], - z=[1, 0], + c=None, + z=None, topboundary="conf", hstar=None, npor=0.3, order=3, ndeg=3, - layers=[0], + layers=None, ): """Element to simulate a building pit with an impermeable wall in ModelMaq. @@ -554,6 +562,12 @@ def __init__( layers: list or np.array layers in which impermeable wall is present. """ + if layers is None: + layers = [0] + if z is None: + z = [1, 0] + if c is None: + c = [] (kaq, c, npor, ltype) = param_maq(kaq, z, c, npor, topboundary) super().__init__( model=model, @@ -577,7 +591,7 @@ def __init__( xy, kaq=1.0, kzoverkh=1.0, - z=[1, 0], + z=None, topboundary="conf", topres=0, topthick=0, @@ -585,7 +599,7 @@ def __init__( npor=0.3, order=3, ndeg=3, - layers=[0], + layers=None, ): """Element to simulate a building pit with an impermeable wall in Model3D. @@ -633,6 +647,10 @@ def __init__( layers: list or np.array layers in which impermeable wall is present. """ + if layers is None: + layers = [0] + if z is None: + z = [1, 0] (kaq, c, npor, ltype) = param_3d(kaq, z, kzoverkh, npor, topboundary, topres) if topboundary == "semi": z = np.hstack((z[0] + topthick, z)) @@ -664,7 +682,7 @@ def __init__( hstar=None, order=3, ndeg=3, - layers=[0], + layers=None, res=np.inf, ): """Element to simulate a building pit with a leaky wall. @@ -713,7 +731,8 @@ def __init__( shape (n_segments,) or (n_layers, n_segments). Default is np.inf, which simulates an impermeable wall. """ - + if layers is None: + layers = [0] super().__init__( model, xy, @@ -744,7 +763,9 @@ def __init__( if np.any(self.res < self.tiny): warn( f"Found resistances smaller than {self.tiny}, " - f"these were replaced by {self.tiny}." + f"these were replaced by {self.tiny}.", + category=UserWarning, + stacklevel=1, ) self.res[self.res < self.tiny] = self.tiny @@ -895,16 +916,22 @@ def __init__( model, xy, kaq=1.0, - z=[1, 0], - c=[], + z=None, + c=None, npor=0.3, topboundary="conf", hstar=None, order=3, ndeg=3, - layers=[0], + layers=None, res=np.inf, ): + if layers is None: + layers = [0] + if c is None: + c = [] + if z is None: + z = [1, 0] (kaq, c, npor, ltype) = param_maq(kaq, z, c, npor, topboundary) super().__init__( model=model, @@ -928,7 +955,7 @@ def __init__( model, xy, kaq=1.0, - z=[1, 0], + z=None, kzoverkh=1.0, npor=0.3, topboundary="conf", @@ -937,7 +964,7 @@ def __init__( hstar=None, order=3, ndeg=3, - layers=[0], + layers=None, res=np.inf, ): """Element to simulate a building pit with a leaky wall in Model3D. @@ -990,6 +1017,10 @@ def __init__( shape (n_segments,) or (n_segments, n_layers). Default is np.inf, which simulates an impermeable wall. """ + if layers is None: + layers = [0] + if z is None: + z = [1, 0] (kaq, c, npor, ltype) = param_3d(kaq, z, kzoverkh, npor, topboundary, topres) if topboundary == "semi": z = np.hstack((z[0] + topthick, z)) diff --git a/timml/inhomogeneity1d.py b/timml/inhomogeneity1d.py index 125892c..48171e6 100644 --- a/timml/inhomogeneity1d.py +++ b/timml/inhomogeneity1d.py @@ -24,7 +24,7 @@ def __init__(self, model, x1, x2, kaq, c, z, npor, ltype, hstar, N): self.addlinesinks = True # Set to False not to add line-sinks def __repr__(self): - return "Inhom1D: " + str(list([self.x1, self.x2])) + return "Inhom1D: " + str([self.x1, self.x2]) def isinside(self, x, y): return (x >= self.x1) and (x < self.x2) @@ -124,13 +124,17 @@ def __init__( x1, x2, kaq=1, - z=[1, 0], - c=[], + z=None, + c=None, npor=0.3, topboundary="conf", hstar=None, N=None, ): + if c is None: + c = [] + if z is None: + z = [1, 0] self.storeinput(inspect.currentframe()) ( kaq, @@ -194,7 +198,7 @@ def __init__( x1, x2, kaq, - z=[1, 0], + z=None, kzoverkh=1, npor=0.3, topboundary="conf", @@ -203,6 +207,8 @@ def __init__( topthick=0.0, N=None, ): + if z is None: + z = [1, 0] self.storeinput(inspect.currentframe()) ( kaq, diff --git a/timml/intlinesink.py b/timml/intlinesink.py index d052276..149a111 100644 --- a/timml/intlinesink.py +++ b/timml/intlinesink.py @@ -262,12 +262,14 @@ def setparams(self, sol): class LeakyIntHeadDiffLineSink(LineSinkHoBase, IntLeakyWallEquation): - """Element to set numerically integrated head along linesink to equal to: + """ + Element to set a numerically integrated head along a linesink. + + The numerically integrated head along the linesink is equal to: Qnormal = H * (headin - headout) / res Used in LeakyBuildingPit element - """ def __init__( diff --git a/timml/linedoublet.py b/timml/linedoublet.py index d14366d..e4e2f49 100644 --- a/timml/linedoublet.py +++ b/timml/linedoublet.py @@ -178,7 +178,6 @@ class ImpLineDoublet(LineDoubletHoBase, DisvecEquation): Parameters ---------- - model : Model object Model to which the element is added x1 : scalar @@ -201,7 +200,6 @@ class ImpLineDoublet(LineDoubletHoBase, DisvecEquation): See Also -------- - :class:`.ImpLineDoubletString` """ @@ -251,7 +249,6 @@ class LeakyLineDoublet(LineDoubletHoBase, LeakyWallEquation): Parameters ---------- - model : Model object Model to which the element is added x1 : scalar @@ -276,7 +273,6 @@ class LeakyLineDoublet(LineDoubletHoBase, LeakyWallEquation): See Also -------- - :class:`.LeakyLineDoubletString` """ @@ -422,7 +418,6 @@ class ImpLineDoubletString(LineDoubletStringBase, DisvecEquation): Parameters ---------- - model : Model object Model to which the element is added xy : array or list @@ -440,11 +435,12 @@ class ImpLineDoubletString(LineDoubletStringBase, DisvecEquation): See Also -------- - :class:`.ImpLineDoublet` """ - def __init__(self, model, xy=[(-1, 0), (1, 0)], layers=0, order=0, label=None): + def __init__(self, model, xy=None, layers=0, order=0, label=None): + if xy is None: + xy = [(-1, 0), (1, 0)] self.storeinput(inspect.currentframe()) LineDoubletStringBase.__init__( self, @@ -473,7 +469,6 @@ class LeakyLineDoubletString(LineDoubletStringBase, LeakyWallEquation): Parameters ---------- - model : Model object Model to which the element is added xy : array or list @@ -493,13 +488,12 @@ class LeakyLineDoubletString(LineDoubletStringBase, LeakyWallEquation): See Also -------- - :class:`.ImpLineDoublet` """ - def __init__( - self, model, xy=[(-1, 0), (1, 0)], res=np.inf, layers=0, order=0, label=None - ): + def __init__(self, model, xy=None, res=np.inf, layers=0, order=0, label=None): + if xy is None: + xy = [(-1, 0), (1, 0)] self.storeinput(inspect.currentframe()) LineDoubletStringBase.__init__( self, diff --git a/timml/linedoublet1d.py b/timml/linedoublet1d.py index c34c70f..3df44b2 100644 --- a/timml/linedoublet1d.py +++ b/timml/linedoublet1d.py @@ -127,7 +127,6 @@ class LeakyLineDoublet1D(LineDoublet1D, LeakyWallEquation): Parameters ---------- - model : Model object Model to which the element is added xld : scalar diff --git a/timml/linesink.py b/timml/linesink.py index b8d1b2c..e1c13c9 100644 --- a/timml/linesink.py +++ b/timml/linesink.py @@ -6,7 +6,7 @@ from . import bessel from .controlpoints import controlpoints, strengthinf_controlpoints from .element import Element -from .equation import HeadEquation, PotentialEquation +from .equation import HeadEquation __all__ = [ "LineSinkBase", @@ -486,7 +486,6 @@ class HeadLineSink(LineSinkHoBase, HeadEquation): Parameters ---------- - model : Model object Model to which the element is added x1 : scalar @@ -520,7 +519,6 @@ class HeadLineSink(LineSinkHoBase, HeadEquation): See Also -------- - :class:`.HeadLineSinkString` """ @@ -589,7 +587,6 @@ class LineSinkDitch(HeadLineSink): Parameters ---------- - model : Model object Model to which the element is added x1 : scalar @@ -620,7 +617,6 @@ class LineSinkDitch(HeadLineSink): See Also -------- - :class:`.LineSinkDitchString` """ @@ -821,9 +817,9 @@ def plot(self, layer=None): class HeadLineSinkStringOLd(LineSinkStringBase, HeadEquation): - def __init__( - self, model, xy=[(-1, 0), (1, 0)], hls=0.0, layers=0, order=0, label=None - ): + def __init__(self, model, xy=None, hls=0.0, layers=0, order=0, label=None): + if xy is None: + xy = [(-1, 0), (1, 0)] self.storeinput(inspect.currentframe()) LineSinkStringBase.__init__( self, @@ -972,7 +968,6 @@ def discharge_per_linesink(self): rv : np.array array of shape (nlay, nlinesinks) """ - Qls = self.parameters[:, 0] * self.dischargeinf() Qls.shape = (self.nls, self.nlayers, self.order + 1) Qls = Qls.sum(axis=2) @@ -983,7 +978,6 @@ def discharge_per_linesink(self): def discharge(self): """Discharge of the element in each layer.""" - rv = np.zeros(self.aq[0].naq) Qls = self.parameters[:, 0] * self.dischargeinf() Qls.shape = (self.nls, self.nlayers, self.order + 1) @@ -1019,7 +1013,6 @@ class HeadLineSinkString(LineSinkStringBase2): Parameters ---------- - model : Model object Model to which the element is added xy : array or list @@ -1051,14 +1044,13 @@ class HeadLineSinkString(LineSinkStringBase2): See Also -------- - :class:`.HeadLineSink` """ def __init__( self, model, - xy=[(-1, 0), (1, 0)], + xy=None, hls=0, res=0, wh=1, @@ -1067,6 +1059,8 @@ def __init__( label=None, name="HeadLineSinkString", ): + if xy is None: + xy = [(-1, 0), (1, 0)] self.storeinput(inspect.currentframe()) LineSinkStringBase2.__init__( self, @@ -1175,7 +1169,6 @@ class LineSinkDitchString(HeadLineSinkString): Parameters ---------- - model : Model object Model to which the element is added xy : array or list @@ -1201,14 +1194,13 @@ class LineSinkDitchString(HeadLineSinkString): See Also -------- - :class:`.LineSinkDitch` """ def __init__( self, model, - xy=[(-1, 0), (1, 0)], + xy=None, Qls=1, res=0, wh=1, @@ -1216,6 +1208,8 @@ def __init__( layers=0, label=None, ): + if xy is None: + xy = [(-1, 0), (1, 0)] self.storeinput(inspect.currentframe()) HeadLineSinkString.__init__( self, @@ -1260,7 +1254,7 @@ class LineSinkContainer(Element): Container class for bunch of line-sinks Required attributes: lslist: list of line-sinks - nls: total number of line-sinks + nls: total number of line-sinks. """ def __init__( @@ -1394,7 +1388,6 @@ class HeadLineSinkContainer(LineSinkContainer): Parameters ---------- - model : Model object Model to which the element is added xydict : dictionary @@ -1421,22 +1414,25 @@ class HeadLineSinkContainer(LineSinkContainer): See Also -------- - :class:`.HeadLineSink` """ def __init__( self, model, - xydict={0: [(-1, 0), (1, 0)]}, + xydict=None, hls=0, res=0, wh=1, order=0, - laydict={0: 0}, + laydict=None, label=None, name="HeadLineSinkContainer", ): + if laydict is None: + laydict = {0: 0} + if xydict is None: + xydict = {0: [(-1, 0), (1, 0)]} self.storeinput(inspect.currentframe()) LineSinkContainer.__init__( self, model, layers=0, order=order, name=name, label=label, aq=None @@ -1481,9 +1477,9 @@ def initialize(self): self.nls = len(self.lslist) for i in range(self.nls): if self.label is not None: - lslabel = self.label + "_" + str(i) + self.label + "_" + str(i) else: - lslabel = self.label + pass LineSinkContainer.initialize(self) def setparams(self, sol): diff --git a/timml/linesink1d.py b/timml/linesink1d.py index 1eb13c0..17670fa 100644 --- a/timml/linesink1d.py +++ b/timml/linesink1d.py @@ -124,7 +124,6 @@ class LineSink1D(LineSink1DBase, MscreenWellEquation): Parameters ---------- - model : Model object Model to which the element is added xls : scalar @@ -179,7 +178,6 @@ class HeadLineSink1D(LineSink1DBase, HeadEquation): Parameters ---------- - model : Model object Model to which the element is added xls : scalar @@ -228,7 +226,7 @@ def setparams(self, sol): class HeadDiffLineSink1D(LineSink1DBase, HeadDiffEquation): - """HeadDiffLineSink1D for left side (xcout)""" + """HeadDiffLineSink1D for left side (xcout).""" def __init__(self, model, xls, label=None, aq=None, aqin=None, aqout=None): LineSink1DBase.__init__( @@ -263,7 +261,7 @@ def setparams(self, sol): class FluxDiffLineSink1D(LineSink1DBase, DisvecDiffEquation): - """HeadDiffLineSink1D for left side (xcout)""" + """HeadDiffLineSink1D for left side (xcout).""" def __init__(self, model, xls, label=None, aq=None, aqin=None, aqout=None): LineSink1DBase.__init__( diff --git a/timml/model.py b/timml/model.py index e7ea023..2e83070 100644 --- a/timml/model.py +++ b/timml/model.py @@ -2,7 +2,6 @@ import inspect # Used for storing the input import multiprocessing as mp -import sys import numpy as np from scipy.integrate import quad_vec @@ -65,7 +64,6 @@ def add_element(self, e): def remove_element(self, e): """Remove element `e` from model.""" - if e.label is not None: self.elementdict.pop(e.label) self.elementlist.remove(e) @@ -86,15 +84,13 @@ def potential(self, x, y, aq=None): return rv def disvec(self, x, y, aq=None): - """Discharge vector at `x`, `y` + """Discharge vector at `x`, `y`. Returns ------- - qxqy : array size (2, naq) first row is Qx in each aquifer layer, second row is Qy """ - if aq is None: aq = self.aq.find_aquifer_data(x, y) rv = np.zeros((2, aq.naq)) @@ -208,7 +204,6 @@ def intnormflux(self, xy, method="legendre", ndeg=10): >>> np.sum(Qn, axis=0) """ - xy = np.array(xy) # convert to array if np.all(xy[-1] == xy[0]): Nsides = len(xy) - 1 @@ -233,16 +228,14 @@ def qztop(self, x, y, aq=None): return rv def head(self, x, y, layers=None, aq=None): - """Head at `x`, `y` + """Head at `x`, `y`. Returns ------- - h : array length `naq` or `len(layers)` head in all `layers` (if not `None`), or all layers of aquifer (otherwise) """ - if aq is None: aq = self.aq.find_aquifer_data(x, y) rv = self.potential(x, y, aq) / aq.T @@ -269,12 +262,10 @@ def headgrid(self, xg, yg, layers=None, printrow=False): ------- h : array size `nlayers, ny, nx` - See also + See Also -------- - :func:`~timml.model.Model.headgrid2` """ - nx, ny = len(xg), len(yg) if layers is None: Nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq @@ -308,12 +299,10 @@ def headgrid2(self, x1, x2, nx, y1, y2, ny, layers=None, printrow=False): ------- h : array size `nlayers, ny, nx` - See also + See Also -------- - :func:`~timml.model.Model.headgrid` """ - xg, yg = np.linspace(x1, x2, nx), np.linspace(y1, y2, ny) return self.headgrid(xg, yg, layers=layers, printrow=printrow) @@ -333,7 +322,6 @@ def headalongline(self, x, y, layers=None): ------- h : array size `nlayers, nx` """ - xg, yg = np.atleast_1d(x), np.atleast_1d(y) if layers is None: Nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq @@ -395,9 +383,12 @@ def velocity(self, x, y, z): def velocomp(self, x, y, z, aq=None, layer_ltype=None): if aq is None: aq = self.aq.find_aquifer_data(x, y) - assert z <= aq.z[0] and z >= aq.z[-1], "z value not inside aquifer" + + if (z > aq.z[0]) or z < (aq.z[-1]): + raise ValueError("z value not inside aquifer") + if layer_ltype is None: - layer, ltype, dummy = aq.findlayer(z) + layer, ltype, _ = aq.findlayer(z) else: layer, ltype = layer_ltype h = self.head(x, y, aq=aq) @@ -616,7 +607,11 @@ class ModelMaq(Model): >>> ml = ModelMaq(kaq=[10, 20], z=[20, 12, 10, 0], c=1000) """ - def __init__(self, kaq=1, z=[1, 0], c=[], npor=0.3, topboundary="conf", hstar=None): + def __init__(self, kaq=1, z=None, c=None, npor=0.3, topboundary="conf", hstar=None): + if c is None: + c = [] + if z is None: + z = [1, 0] self.storeinput(inspect.currentframe()) kaq, c, npor, ltype = param_maq(kaq, z, c, npor, topboundary) Model.__init__(self, kaq, c, z, npor, ltype) @@ -670,7 +665,7 @@ class Model3D(Model): def __init__( self, kaq=1, - z=[1, 0], + z=None, kzoverkh=1, npor=0.3, topboundary="conf", @@ -682,7 +677,10 @@ def __init__( for semi-confined aquifers, set top equal to 'semi' and provide topres: resistance of top tophick: thickness of top - hstar: head above top""" + hstar: head above top. + """ + if z is None: + z = [1, 0] self.storeinput(inspect.currentframe()) kaq, c, npor, ltype = param_3d(kaq, z, kzoverkh, npor, topboundary, topres) if topboundary == "semi": diff --git a/timml/stripareasink.py b/timml/stripareasink.py index 20ec7f2..48b0129 100644 --- a/timml/stripareasink.py +++ b/timml/stripareasink.py @@ -1,5 +1,3 @@ -import inspect # Used for storing the input - import numpy as np from .element import Element @@ -120,9 +118,6 @@ def changetrace( u = u1 * (1.0 + eps) # Go just beyond circle else: u = u2 * (1.0 + eps) # Go just beyond circle - xn = x1 + u * (x2 - x1) - yn = y1 + u * (y2 - y1) - zn = xyzt1[2] + u * (xyzt2[2] - xyzt1[2]) xyztnew = xyzt1 + u * (xyzt2 - xyzt1) return changed, terminate, xyztnew, message @@ -250,8 +245,6 @@ def changetrace( u = u1 * (1.0 + eps) # Go just beyond circle else: u = u2 * (1.0 + eps) # Go just beyond circle - xn = x1 + u * (x2 - x1) - yn = y1 + u * (y2 - y1) - zn = xyzt1[2] + u * (xyzt2[2] - xyzt1[2]) + xyzt1[2] + u * (xyzt2[2] - xyzt1[2]) xyztnew = xyzt1 + u * (xyzt2 - xyzt1) return changed, terminate, xyztnew, message diff --git a/timml/trace.py b/timml/trace.py index 2da4141..647a974 100644 --- a/timml/trace.py +++ b/timml/trace.py @@ -20,7 +20,7 @@ def timtraceline( vstepfrac=0.2, tmax=1e12, nstepmax=100, - win=[-1e30, 1e30, -1e30, 1e30], + win=None, silent=False, returnlayers=False, *, @@ -59,6 +59,8 @@ def timtraceline( - "message": termination message - "complete": True if terminated correctly """ + if win is None: + win = [-1e30, 1e30, -1e30, 1e30] verbose = False # used for debugging if not metadata: warnings.warn(_future_warning_metadata, FutureWarning, stacklevel=2) @@ -288,7 +290,7 @@ def timtracelines( tmax=1e12, nstepmax=100, silent=".", - win=[-1e30, 1e30, -1e30, 1e30], + win=None, *, metadata=False, ): @@ -320,6 +322,8 @@ def timtracelines( if False, return list of xyzt arrays if True, return list of result dicionaries """ + if win is None: + win = [-1e30, 1e30, -1e30, 1e30] xyztlist = [] for x, y, z in zip(xstart, ystart, zstart): xyztlist.append( diff --git a/timml/util.py b/timml/util.py index c58465c..9e7328f 100644 --- a/timml/util.py +++ b/timml/util.py @@ -2,7 +2,7 @@ import numpy as np from matplotlib.collections import LineCollection -from .trace import timtraceline, timtracelines +from .trace import timtraceline plt.rcParams["contour.negative_linestyle"] = "solid" @@ -24,7 +24,6 @@ def plot( Parameters ---------- - win : list or tuple [xmin, xmax, ymin, ymax] newfig : boolean (default True) @@ -41,10 +40,8 @@ def plot( Returns ------- - None """ - if newfig: plt.figure(figsize=figsize) ax1 = None @@ -120,7 +117,6 @@ def contour( Parameters ---------- - win : list or tuple [xmin, xmax, ymin, ymax] ngr : scalar, tuple or list @@ -149,10 +145,8 @@ def contour( Returns ------- - cs : list of contour sets for each contoured layer """ - x1, x2, y1, y2 = win if np.isscalar(ngr): nx = ny = ngr @@ -167,9 +161,9 @@ def contour( # color if color is None: c = plt.rcParams["axes.prop_cycle"].by_key()["color"] - elif type(color) is str: + elif isinstance(color, str): c = len(layers) * [color] - elif type(color) is list: + elif isinstance(color, list): c = color if len(c) < len(layers): n = np.ceil(self.aq.naq / len(c)) @@ -185,7 +179,7 @@ def contour( if labels: fmt = "%1." + str(decimals) + "f" plt.clabel(cs, fmt=fmt) - if type(legend) is list: + if isinstance(legend, list): plt.legend(cshandlelist, legend) elif legend: legendlist = ["layer " + str(i) for i in layers] @@ -213,7 +207,6 @@ def vcontour( Parameters ---------- - win : list or tuple [xmin, xmax, ymin, ymax] n : integer @@ -240,10 +233,8 @@ def vcontour( Returns ------- - cs : contour set """ - x1, x2, y1, y2 = win h = self.headalongline( np.linspace(x1 + nudge, x2 - nudge, n), @@ -283,7 +274,7 @@ def tracelines( silent=".", color=None, orientation="hor", - win=[-1e30, 1e30, -1e30, 1e30], + win=None, newfig=False, figsize=None, *, @@ -329,15 +320,16 @@ def tracelines( Returns ------- - traces : result only if return_traces = True """ + if win is None: + win = [-1e30, 1e30, -1e30, 1e30] if color is None: c = plt.rcParams["axes.prop_cycle"].by_key()["color"] - elif type(color) is str: + elif isinstance(color, str): c = self.aq.naq * [color] - elif type(color) is list: + elif isinstance(color, list): c = color if len(c) < self.aq.naq: n = int(np.ceil(self.aq.naq / len(c))) @@ -427,7 +419,6 @@ def vcontoursf1D( Parameters ---------- - x1 : scalar left edge of contour domain x2 : scalar @@ -455,7 +446,6 @@ def vcontoursf1D( Returns ------- - ax : axis """ naq = self.aq.naq diff --git a/timml/well.py b/timml/well.py index 1c05107..6cfb4e4 100644 --- a/timml/well.py +++ b/timml/well.py @@ -5,7 +5,7 @@ from scipy.special import k0, k1 from .element import Element -from .equation import MscreenWellEquation, PotentialEquation, MscreenWellNoflowEquation +from .equation import MscreenWellEquation, MscreenWellNoflowEquation, PotentialEquation from .trace import timtracelines __all__ = ["WellBase", "Well", "HeadWell"] @@ -115,7 +115,6 @@ def headinside(self): array (length number of screens) Head inside the well for each screen """ - h = self.model.head(self.xw + self.rw, self.yw, layers=self.layers) return h - self.resfac * self.parameters[:, 0] @@ -127,7 +126,6 @@ def discharge(self): array (length number of layers) Discharge in each screen with zeros for layers that are not screened """ - Q = np.zeros(self.aq.naq) Q[self.layers] = self.parameters[:, 0] return Q @@ -240,7 +238,7 @@ def plotcapzone( silent=".", color=None, orientation="hor", - win=[-1e30, 1e30, -1e30, 1e30], + win=None, newfig=False, figsize=None, *, @@ -276,6 +274,8 @@ def plotcapzone( figsize : tuple of integers, optional, default: None width, height in inches. """ + if win is None: + win = [-1e30, 1e30, -1e30, 1e30] if not return_traces: metadata = True # suppress future warning from timtraceline xstart, ystart, zstart = self.capzonestart(nt, zstart) @@ -301,7 +301,7 @@ def plotcapzone( class Well(WellBase, MscreenWellEquation): - """Well Class to create a well with a specified discharge. + r"""Well Class to create a well with a specified discharge. Notes ----- @@ -388,7 +388,7 @@ def setparams(self, sol): class HeadWell(WellBase, PotentialEquation): - """HeadWell Class to create a well with a specified head inside the well. + r"""HeadWell Class to create a well with a specified head inside the well. Notes ----- @@ -586,5 +586,5 @@ def disvecinf(self, x, y, aq=None): qy[:] = -kone * yminyw / (r * aq.lab) / (2 * np.pi) rv[0] = self.aq.coef[self.layers] * qx rv[1] = self.aq.coef[self.layers] * qy - #rv[0], rv[1] = qx, qy + # rv[0], rv[1] = qx, qy return rv