diff --git a/.docs/conf.py b/.docs/conf.py index 4639f4dd..54ea5cac 100644 --- a/.docs/conf.py +++ b/.docs/conf.py @@ -20,52 +20,27 @@ sys.path.insert(0, os.path.abspath("..")) from flopy import __version__ import pymake +# -- determine if running on readthedocs ------------------------------------ +on_rtd = os.environ.get('READTHEDOCS') == 'True' + # -- get the MODFLOW executables -------------------------------------------- -ws = ".bin" -if os.path.isdir(ws): - shutil.rmtree(ws) -os.makedirs(ws) -osname = sys.platform.lower() -if osname == "darwin": - platform = "mac" -else: - platform = osname -pymake.getmfexes(pth=ws, platform=platform, verbose=True) +if not on_rtd: + ws = ".bin" + if os.path.isdir(ws): + shutil.rmtree(ws) + os.makedirs(ws) + osname = sys.platform.lower() + if osname == "darwin": + platform = "mac" + else: + platform = osname + pymake.getmfexes(pth=ws, platform=platform, verbose=True) -# -- copy the example problems ---------------------------------------------- -srcdir = "pysrc" -ws = ".working" -if os.path.isdir(ws): - shutil.rmtree(ws) -os.makedirs(ws) -scripts = [] -for fname in os.listdir(srcdir): - if fname.endswith(".py"): - scripts.append(fname) - src = os.path.join(srcdir, fname) - dst = os.path.join(ws, fname) - shutil.copyfile(src, dst) - -# -- run the example problems ----------------------------------------------- -for s in scripts: - args = ("python", s) - proc = Popen(args, stdout=PIPE, stderr=PIPE, cwd=ws) - stdout, stderr = proc.communicate() - if stdout: - print(stdout.decode("utf-8")) - if stderr: - print("Errors:\n{}".format(stderr.decode("utf-8"))) - -# -- copy the output to _static directory ----------------------------------- -ws_dst = "_static" -if os.path.isdir(ws_dst): - shutil.rmtree(ws_dst) -os.makedirs(ws_dst) -for fname in os.listdir(ws): - if fname.endswith(".png"): - src = os.path.join(ws, fname) - dst = os.path.join(ws_dst, fname) - shutil.copyfile(src, dst) +# -- convert the tutorial scripts ------------------------------------------- +if not on_rtd: + cmd = ("python", "tutorials2ipynb.py") + print(" ".join(cmd)) + os.system(" ".join(cmd)) # -- Create the flopy rst files --------------------------------------------- args = ("sphinx-apidoc", "-e", "-o", "source/", "../flopy/") @@ -78,7 +53,6 @@ if stderr: # -- Project information ----------------------------------------------------- - project = "flopy" copyright = "2020, Bakker, Mark, Post, Vincent, Langevin, C. D., Hughes, J. D., White, J. T., Leaf, A. T., Paulinski, S. R., Larsen, J. D., Toews, M. W., Morway, E. D., Bellino, J. C., Starn, J. J., and Fienen, M. N." author = "Bakker, Mark, Post, Vincent, Langevin, C. D., Hughes, J. D., White, J. T., Leaf, A. T., Paulinski, S. R., Larsen, J. D., Toews, M. W., Morway, E. D., Bellino, J. C., Starn, J. J., and Fienen, M. N." @@ -111,6 +85,18 @@ extensions = [ "recommonmark", ] +# nbsphinx_execute_arguments = [ +# "--InlineBackend.figure_formats={'svg', 'pdf'}", +# "--InlineBackend.rc={'figure.dpi': 200}", +# ] + +# Settings for GitHub actions integration +if on_rtd: + extensions.append("rtds_action") + rtds_action_github_repo = "modflowpy/flopy" + rtds_action_path = "_notebooks" + rtds_action_artifact_prefix = "notebooks-for-" + rtds_action_github_token = os.environ.get("GITHUB_TOKEN", None) # Add any paths that contain templates here, relative to this directory. templates_path = ["_templates"] @@ -183,7 +169,7 @@ html_css_files = [ # A shorter title for the navigation bar. Default is the same as html_title. html_short_title = "flopy" -html_favicon = ".._images/flopylogo.png" +html_favicon = "_images/flopylogo.png" # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, diff --git a/.docs/local_test.sh b/.docs/local_test.sh deleted file mode 100644 index dcb06860..00000000 --- a/.docs/local_test.sh +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash - -sphinx-apidoc -e -o source/ ../flopy/ -make html diff --git a/.docs/tutorial1.rst b/.docs/tutorial1.rst deleted file mode 100644 index 12d6f02e..00000000 --- a/.docs/tutorial1.rst +++ /dev/null @@ -1,202 +0,0 @@ -Tutorial 1: Confined Steady-State Flow Model -============================================ - -This tutorial demonstrates use of FloPy to develop a simple MODFLOW model. -Note that you can access the latest version this tutorial python script from -`here `_. - -Getting Started ---------------- -If FloPy has been properly installed, then it can be imported as follows:: - - - import flopy - - -Now that we can import flopy, we begin creating our simple MODFLOW model. - -Creating the MODFLOW Model --------------------------- -One of the nice things about creating models in python is that it is very easy -to change one or two things and completely change the grid resolution for your -model. So in this example, we will design our python script so that the number -of layers, columns, and rows can be easily changed. - -We can create a very simple MODFLOW model that has a basic package (BAS), -discretization input file (DIS), layer-property flow (LPF) package, output -control (OC), and preconditioned conjugate gradient (PCG) solver. Each one of -these has its own input file, which will be created automatically by flopy, -provided that we pass flopy the correct information. - -Discretization -^^^^^^^^^^^^^^ - -We start by creating our flopy model object as follows:: - - - # Assign name and create modflow model object - modelname = 'tutorial1' - mf = flopy.modflow.Modflow(modelname, exe_name='mf2005') - -Next, let's proceed by defining our model domain and creating a MODFLOW grid -to span the domain:: - - # Model domain and grid definition - Lx = 1000. - Ly = 1000. - ztop = 0. - zbot = -50. - nlay = 1 - nrow = 10 - ncol = 10 - delr = Lx / ncol - delc = Ly / nrow - delv = (ztop - zbot) / nlay - botm = np.linspace(ztop, zbot, nlay + 1) - -With this information, we can now create the flopy discretization object by -entering the following:: - - # Create the discretization object - dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, - top=ztop, botm=botm[1:]) - -The obvious question at this point is, how do I know which arguments are -required by this strange thing called flopy.modflow.ModflowDis? Fortunately, -there is an online help page for each one of the model objects. The page for -the DIS input file is located at `flopy.modflow.mfdis `__. - -Basic Package -^^^^^^^^^^^^^ - -Next we can create a flopy object that represents the MODFLOW Basic Package. -Details on the flopy BAS class are at: `flopy.modflow.mfbas `__. -For this simple model, we will assign constant head values of 10. and 0. to the -first and last model columns (in all layers), respectively. The python code -for doing this is:: - - # Variables for the BAS package - ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) - ibound[:, :, 0] = -1 - ibound[:, :, -1] = -1 - strt = np.ones((nlay, nrow, ncol), dtype=np.float32) - strt[:, :, 0] = 10. - strt[:, :, -1] = 0. - bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) - -Layer-Property Flow Package -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Details on the flopy LPF class are at: `flopy.modflow.mflpf `__. -Values of 10. are assigned for the horizontal and vertical hydraulic -conductivity:: - - # Add LPF package to the MODFLOW model - lpf = flopy.modflow.ModflowLpf(mf, hk=10., vka=10., ipakcb=53) - -Because we did not specify a value for laytyp, Flopy will use the default value -of 0, which means that this model will be confined. - -Output Control -^^^^^^^^^^^^^^ - -Details on the flopy OC class are at: `flopy.modflow.mfoc `__. Here -we can use the default OC settings by specifying the following:: - - # Add OC package to the MODFLOW model - spd = {(0, 0): ['print head', 'print budget', 'save head', 'save budget']} - oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True) - -The stress period dictionary is used to set what output is saved for the -corresponding stress period and time step. In this case, the tuple (0, 0) -means that stress period 1 and time step 1 for MODFLOW will have output saved. -Head and budgets will be printed and head and budget information will be saved. - -Preconditioned Conjugate Gradient Package -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Details on the flopy PCG class are at: `flopy.modflow.mfpcg `__. -The default settings used by flopy will be used by specifying the following -commands:: - - # Add PCG package to the MODFLOW model - pcg = flopy.modflow.ModflowPcg(mf) - -Writing the MODFLOW Data Files -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -The MODFLOW input data files are written by simply issuing the following:: - - # Write the MODFLOW model input files - mf.write_input() - -Running the Modeling --------------------- - -Flopy can also be used to run the model. The model object (mf in this example) -has an attached method that will run the model. For this to work, the MODFLOW -program must be located somewhere within the system path, or within the working -directory. In this example, we have specified that the name of the executable -program is 'mf2005'. Issue the following to run the model:: - - # Run the MODFLOW model - success, buff = mf.run_model() - -Here we have used run_model, and we could also have specified values for the -optional keywords silent, pause, and report. - -Post-Processing the Results ---------------------------- - -Now that we have successfully built and run our MODFLOW model, we can look at -the results. MODFLOW writes the simulated heads to a binary data output file. -We cannot look at these heads with a text editor, but flopy has a binary -utility that can be used to read the heads. The following statements will -read the binary head file and create a plot of simulated heads for layer 1:: - - import matplotlib.pyplot as plt - import flopy.utils.binaryfile as bf - plt.subplot(1, 1, 1, aspect='equal') - hds = bf.HeadFile(modelname + '.hds') - head = hds.get_data(totim=1.0) - levels = np.arange(1, 10, 1) - extent = (delr / 2., Lx - delr / 2., Ly - delc / 2., delc / 2.) - plt.contour(head[0, :, :], levels=levels, extent=extent) - plt.savefig('tutorial1a.png') - -If everything has worked properly, you should see the following head contours. - -.. figure:: _static/tutorial1a.png - :alt: head contours in first layer - :scale: 100 % - :align: left - -Flopy also has some pre-canned plotting capabilities can can be accessed using -the PlotMapView class. The following code shows how to use the plotmapview -class to plot boundary conditions (IBOUND), plot the grid, plot head contours, -and plot vectors:: - - fig = plt.figure(figsize=(10,10)) - ax = fig.add_subplot(1, 1, 1, aspect='equal') - - hds = bf.HeadFile(modelname + '.hds') - times = hds.get_times() - head = hds.get_data(totim=times[-1]) - levels = np.linspace(0, 10, 11) - - cbb = bf.CellBudgetFile(modelname + '.cbc') - kstpkper_list = cbb.get_kstpkper() - frf = cbb.get_data(text='FLOW RIGHT FACE', totim=times[-1])[0] - fff = cbb.get_data(text='FLOW FRONT FACE', totim=times[-1])[0] - - pmv = flopy.plot.PlotMapView(model=mf, layer=0) - qm = pmv.plot_ibound() - lc = pmv.plot_grid() - cs = pmv.contour_array(head, levels=levels) - quiver = pmv.plot_discharge(frf, fff, head=head) - plt.savefig('tutorial1b.png') - -.. figure:: _static/tutorial1b.png - :alt: head contours in first layer - :scale: 100 % - :align: left diff --git a/.docs/tutorial2.rst b/.docs/tutorial2.rst deleted file mode 100644 index 5959568b..00000000 --- a/.docs/tutorial2.rst +++ /dev/null @@ -1,295 +0,0 @@ -Tutorial 2: Unconfined Transient Flow Model -=========================================== - -In this example, we will convert the tutorial 1 model into an unconfined, -transient flow model with time varying boundaries. Instead of using constant -heads for the left and right boundaries (by setting ibound to -1), we will use -general head boundaries. We will have the model consider the following -conditions: - -* Initial conditions -- head is 10.0 everywhere -* Period 1 (1 day) -- steady state with left and right GHB stage = 10. -* Period 2 (100 days) -- left GHB with stage = 10., right GHB with stage set - to 0. -* Period 3 (100 days) -- pumping well at model center with rate = -500., left - and right GHB = 10., and 0. - -We will start with selected model commands from the previous tutorial. - -Note that you can access the latest version of this tutorial python script from -`here `_. - -Getting Started ---------------- -As shown in the previous tutorial, import flopy using your preferred method, -such as:: - - import numpy as np - import flopy - -Creating the MODFLOW Model --------------------------- - -Define the Model Extent, Grid Resolution, and Characteristics -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Assign the model information:: - - # Model domain and grid definition - Lx = 1000. - Ly = 1000. - ztop = 10. - zbot = -50. - nlay = 1 - nrow = 10 - ncol = 10 - delr = Lx / ncol - delc = Ly / nrow - delv = (ztop - zbot) / nlay - botm = np.linspace(ztop, zbot, nlay + 1) - hk = 1. - vka = 1. - sy = 0.1 - ss = 1.e-4 - laytyp = 1 - - # Variables for the BAS package - # Note that changes from the previous tutorial! - ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) - strt = 10. * np.ones((nlay, nrow, ncol), dtype=np.float32) - - -Define the Stress Periods -^^^^^^^^^^^^^^^^^^^^^^^^^ - -To create a model with multiple stress periods, we need to define nper, perlen, -nstp, and steady. This is done in the following block in a manner that allows -us to pass these variable directly to the discretization object:: - - # Time step parameters - nper = 3 - perlen = [1, 100, 100] - nstp = [1, 100, 100] - steady = [True, False, False] - -Create Time-Invariant Flopy Objects -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -With this information, we can now create the static flopy objects that do -not change with time:: - - # Flopy objects - modelname = 'tutorial2' - mf = flopy.modflow.Modflow(modelname, exe_name='mf2005') - dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, - top=ztop, botm=botm[1:], - nper=nper, perlen=perlen, nstp=nstp, steady=steady) - bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) - lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, ipakcb=53) - pcg = flopy.modflow.ModflowPcg(mf) - - -Transient General-Head Boundary Package -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -At this point, our model is ready to add our transient boundary packages. -First, we will create the GHB object, which is of the following -type: `flopy.modflow.ModflowGhb `__. - -The key to creating Flopy transient boundary packages is recognizing that -the boundary data is stored in a dictionary with key values equal to the -zero-based stress period number and values equal to the boundary conditions -for that stress period. For a GHB the values can be a two-dimensional nested -list of [layer, row, column, stage, conductance]:: - - # Make list for stress period 1 - stageleft = 10. - stageright = 10. - bound_sp1 = [] - for il in range(nlay): - condleft = hk * (stageleft - zbot) * delc - condright = hk * (stageright - zbot) * delc - for ir in range(nrow): - bound_sp1.append([il, ir, 0, stageleft, condleft]) - bound_sp1.append([il, ir, ncol - 1, stageright, condright]) - print('Adding ', len(bound_sp1), 'GHBs for stress period 1.') - - # Make list for stress period 2 - stageleft = 10. - stageright = 0. - condleft = hk * (stageleft - zbot) * delc - condright = hk * (stageright - zbot) * delc - bound_sp2 = [] - for il in range(nlay): - for ir in range(nrow): - bound_sp2.append([il, ir, 0, stageleft, condleft]) - bound_sp2.append([il, ir, ncol - 1, stageright, condright]) - print('Adding ', len(bound_sp2), 'GHBs for stress period 2.') - - # We do not need to add a dictionary entry for stress period 3. - # Flopy will automatically take the list from stress period 2 and apply it - # to the end of the simulation, if necessary - stress_period_data = {0: bound_sp1, 1: bound_sp2} - - # Create the flopy ghb object - ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=stress_period_data) - - -Transient Well Package -^^^^^^^^^^^^^^^^^^^^^^ - -Now we can create the well package object, which is of the type, -`flopy.modflow.ModflowWel `__.:: - - # Create the well package - # Remember to use zero-based layer, row, column indices! - pumping_rate = -500. - wel_sp1 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]] - wel_sp2 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]] - wel_sp3 = [[0, nrow/2 - 1, ncol/2 - 1, pumping_rate]] - stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3} - wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data) - - -Output Control -^^^^^^^^^^^^^^ - -Here we create the output control package object, which is of the -type `flopy.modflow.ModflowOc `__.:: - - stress_period_data = {} - for kper in range(nper): - for kstp in range(nstp[kper]): - stress_period_data[(kper, kstp)] = ['save head', - 'save drawdown', - 'save budget', - 'print head', - 'print budget'] - oc = flopy.modflow.ModflowOc(mf, stress_period_data=stress_period_data, - compact=True) - - -Running the Modeling --------------------- - -Run the model with the run_model method, which returns a success flag and -the stream of output. With run_model, we have some finer control, that allows us to suppress the output.:: - - # Write the model input files - mf.write_input() - - # Run the model - success, mfoutput = mf.run_model(silent=True, pause=False, report=True) - if not success: - raise Exception('MODFLOW did not terminate normally.') - -Post-Processing the Results ---------------------------- - -Once again, we can read heads from the MODFLOW binary output file, using -the `flopy.utils.binaryfile `__ module. Included with the -HeadFile object are several methods that we will use here: - -* get_times() will return a list of times contained in the binary head file -* get_data() will return a three-dimensional head array for the specified time -* get_ts() will return a time series array [ntimes, headval] for the specified cell - -Using these methods, we can create head plots and hydrographs from the -model results.:: - - # Imports - import matplotlib.pyplot as plt - import flopy.utils.binaryfile as bf - - # Create the headfile and budget file objects - headobj = bf.HeadFile(modelname + '.hds') - times = headobj.get_times() - cbb = bf.CellBudgetFile(modelname + '.cbc') - - # Setup contour parameters - levels = np.linspace(0, 10, 11) - extent = (delr / 2., Lx - delr / 2., delc / 2., Ly - delc / 2.) - print('Levels: ', levels) - print('Extent: ', extent) - - # Well point - wpt = ((float(ncol / 2) - 0.5) * delr, (float(nrow / 2 - 1) + 0.5) * delc) - wpt = (450., 550.) - - # Make the plots - mytimes = [1.0, 101.0, 201.0] - for iplot, time in enumerate(mytimes): - print('*****Processing time: ', time) - head = headobj.get_data(totim=time) - #Print statistics - print('Head statistics') - print(' min: ', head.min()) - print(' max: ', head.max()) - print(' std: ', head.std()) - - # Extract flow right face and flow front face - frf = cbb.get_data(text='FLOW RIGHT FACE', totim=time)[0] - fff = cbb.get_data(text='FLOW FRONT FACE', totim=time)[0] - - #Create the plot - f = plt.figure() - plt.subplot(1, 1, 1, aspect='equal') - plt.title('stress period ' + str(iplot + 1)) - - - pmv = flopy.plot.PlotMapView(model=mf, layer=0) - qm = pmv.plot_ibound() - lc = pmv.plot_grid() - qm = pmv.plot_bc('GHB', alpha=0.5) - cs = pmv.contour_array(head, levels=levels) - plt.clabel(cs, inline=1, fontsize=10, fmt='%1.1f') - quiver = pmv.plot_discharge(frf, fff, head=head) - - - mfc = 'None' - if (iplot + 1) == len(mytimes): - mfc='black' - plt.plot(wpt[0], wpt[1], lw=0, marker='o', markersize=8, - markeredgewidth=0.5, - markeredgecolor='black', markerfacecolor=mfc, zorder=9) - plt.text(wpt[0] + 25, wpt[1] - 25, 'well', size=12, zorder=12) - plt.savefig('tutorial2-{}.png'.format(iplot)) - -If everything has worked properly, you should see the following head contours. - -.. figure:: _static/tutorial2-0.png - :alt: head contours for stress period 1 - :scale: 100 % - :align: center - -.. figure:: _static/tutorial2-1.png - :alt: head contours for stress period 2 - :scale: 100 % - :align: center - -.. figure:: _static/tutorial2-2.png - :alt: head contours for stress period 3 - :scale: 100 % - :align: center - - -Plot Head Versus Time -^^^^^^^^^^^^^^^^^^^^^ -Make a plot of head versus time by extracting the binary heads from -the headobj:: - - # Plot the head versus time - idx = (0, int(nrow / 2) - 1, int(ncol / 2) - 1) - ts = headobj.get_ts(idx) - plt.subplot(1, 1, 1) - ttl = 'Head at cell ({0},{1},{2})'.format(idx[0] + 1, idx[1] + 1, idx[2] + 1) - plt.title(ttl) - plt.xlabel('time') - plt.ylabel('head') - plt.plot(ts[:, 0], ts[:, 1], 'bo-') - plt.savefig('tutorial2-ts.png') - -.. figure:: _static/tutorial2-ts.png - :alt: head contours in first layer - :scale: 100 % - :align: center diff --git a/.docs/tutorials.rst b/.docs/tutorials.rst index 90866fe7..2f89ffd6 100644 --- a/.docs/tutorials.rst +++ b/.docs/tutorials.rst @@ -1,10 +1,33 @@ Tutorials ========= +MODFLOW 6 Tutorials +----------------- + Contents: .. toctree:: :maxdepth: 2 - - tutorial1 - tutorial2 + + _notebooks/tutorial01_mf6 + +MODFLOW Tutorials +----------------- + +Contents: + +.. toctree:: + :maxdepth: 2 + + _notebooks/tutorial01_mf + _notebooks/tutorial02_mf + +SEAWAT Tutorials +----------------- + +Contents: + +.. toctree:: + :maxdepth: 2 + + _notebooks/tutorial01_seawat diff --git a/.docs/tutorials2ipynb.py b/.docs/tutorials2ipynb.py new file mode 100644 index 00000000..91c36dee --- /dev/null +++ b/.docs/tutorials2ipynb.py @@ -0,0 +1,36 @@ +import os +import shutil + +wpth = ".working" +if os.path.isdir(wpth): + shutil.rmtree(wpth) +os.makedirs(wpth) +# copy the python files +pth = os.path.join("..", "examples", "Tutorials") +py_files = [file_name for file_name in os.listdir(pth) + if file_name.endswith(".py")] +for file_name in py_files: + src = os.path.join(pth, file_name) + dst = os.path.join(wpth, file_name) + shutil.copyfile(src, dst) +py_pth = os.path.join(wpth, "*.py") +cmd = ( + "jupytext", + "--to ipynb", + "--execute", + py_pth, + ) +print(" ".join(cmd)) +os.system(" ".join(cmd)) + +npth = "_notebooks" +# copy notebooks +if os.path.isdir(npth): + shutil.rmtree(npth) +os.makedirs(npth) +for file_name in py_files: + src = os.path.join(wpth, file_name.replace(".py", ".ipynb")) + dst = os.path.join(npth, file_name.replace(".py", ".ipynb")) + shutil.copyfile(src, dst) +shutil.rmtree(".working") + diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e353a626..03873d2d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -104,12 +104,18 @@ jobs: activate-environment: flopy use-only-tar-bz2: true - - name: Add jupyter and nbconvert to one matrix entry + - name: Add jupyter and nbconvert to notebooks run if: matrix.run-type == 'nb' shell: bash -l {0} run: | conda install -n flopy -c conda-forge jupyter nbconvert + - name: Add jupyter and jupytext to scripts run + if: matrix.run-type == 'script' + shell: bash -l {0} + run: | + conda install -n flopy -c conda-forge jupyter jupytext + - name: Add packages to flopy environment shell: bash -l {0} run: | @@ -154,3 +160,44 @@ jobs: uses: codecov/codecov-action@v1.0.12 with: file: ./coverage.xml + + - name: Run jupytext on tutorials + if: matrix.run-type == 'script' + shell: bash -l {0} + run: | + cd .docs/ + python tutorials2ipynb.py + cd ../ + + - name: Upload completed jupyter notebooks as an artifact for ReadtheDocs + if: matrix.run-type == 'script' && github.repository_owner == 'modflowpy' + uses: actions/upload-artifact@v2 + with: + name: notebooks-for-${{ github.sha }} + path: | + .docs/_notebooks + + # trigger rtd if previous job was successful + rtd: + name: ReadtheDocs + needs: flopyCI + runs-on: ubuntu-latest + + if: github.repository_owner == 'modflowpy' + steps: + - name: Checkout flopy repo + uses: actions/checkout@v2 + + - name: Output repo information + run: | + echo ${{ github.repository_owner }} + echo ${{ github.repository }} + echo ${{ github.ref }} + + - name: Trigger RTDs build on master and develop branches + if: github.ref == 'refs/heads/master' || github.ref == 'refs/heads/develop' + uses: dfm/rtds-action@v1.0.0 + with: + webhook_url: ${{ secrets.RTDS_WEBHOOK_URL }} + webhook_token: ${{ secrets.RTDS_WEBHOOK_TOKEN }} + commit_ref: ${{ github.ref }} diff --git a/.gitignore b/.gitignore index 9cb0e5bb..44811659 100644 --- a/.gitignore +++ b/.gitignore @@ -85,4 +85,6 @@ autotest/.noseids .docs/_templates .docs/source .docs/.bin +.docs/_temp +.docs/_notebooks diff --git a/autotest/autotest_scripts.py b/autotest/autotest_scripts.py index 5b352082..5ffcc153 100644 --- a/autotest/autotest_scripts.py +++ b/autotest/autotest_scripts.py @@ -3,6 +3,7 @@ from __future__ import print_function import os import sys import shutil +from subprocess import Popen, PIPE exclude = ['flopy_swi2_ex2.py', 'flopy_swi2_ex5.py'] for arg in sys.argv: @@ -10,6 +11,7 @@ for arg in sys.argv: exclude = [] sdir = os.path.join('..', 'examples', 'scripts') +tdir = os.path.join("..", "examples", "Tutorials") # make working directories tempdir = os.path.join('.', 'temp') @@ -17,17 +19,20 @@ if os.path.isdir(tempdir): shutil.rmtree(tempdir) os.mkdir(tempdir) -testdir = os.path.join('.', 'temp', 'scripts') -if os.path.isdir(testdir): - shutil.rmtree(testdir) -os.mkdir(testdir) +testdirs = (os.path.join('.', 'temp', 'scripts'), + os.path.join('.', 'temp', 'tutorials'), + ) +for testdir in testdirs: + if os.path.isdir(testdir): + shutil.rmtree(testdir) + os.mkdir(testdir) -# add testdir to python path -sys.path.append(testdir) + # add testdir to python path + sys.path.append(testdir) -def copy_scripts(): - files = [f for f in os.listdir(sdir) if f.endswith('.py')] +def copy_scripts(src_dir, dst_dir): + files = [f for f in sorted(os.listdir(src_dir)) if f.endswith('.py')] # exclude unwanted files for e in exclude: @@ -35,13 +40,13 @@ def copy_scripts(): files.remove(e) # copy files - for fn in files: - pth = os.path.join(sdir, fn) - opth = os.path.join(testdir, fn) + for file_name in files: + src = os.path.join(src_dir, file_name) + dst = os.path.join(dst_dir, file_name) # copy script - print('copying {} from {} to {}'.format(fn, sdir, testdir)) - shutil.copyfile(pth, opth) + print('copying {} from {} to {}'.format(file_name, src_dir, testdir)) + shutil.copyfile(src, dst) return files @@ -52,7 +57,7 @@ def import_from(mod, name): return main -def run_scripts(fn): +def run_scripts(fn, testdir): # import run function from scripts s = os.path.splitext(fn)[0] run = import_from(s, 'run') @@ -73,19 +78,42 @@ def run_scripts(fn): assert ival == 0, 'could not run {}'.format(fn) -def test_notebooks(): +def run_tutorial_scripts(fn, testdir): + args = ("python", fn) + print("running...'{}'".format(" ".join(args))) + proc = Popen(args, stdout=PIPE, stderr=PIPE, cwd=testdir) + stdout, stderr = proc.communicate() + if stdout: + print(stdout.decode("utf-8")) + if stderr: + print("Errors:\n{}".format(stderr.decode("utf-8"))) + return + + +def test_scripts(): # get list of scripts to run - files = copy_scripts() + files = copy_scripts(sdir, testdirs[0]) for fn in files: - yield run_scripts, fn + yield run_scripts, fn, testdirs[0] + + +def test_tutorial_scripts(): + # get list of scripts to run + files = copy_scripts(tdir, testdirs[1]) + + for fn in files: + yield run_tutorial_scripts, fn, testdirs[1] if __name__ == '__main__': + # # get list of scripts to run + # files = copy_scripts(sdir, testdirs[0]) + # for fn in files: + # run_scripts(fn, testdirs[0]) - # get list of scripts to run - files = copy_scripts() - + # get list of tutorial scripts to run + files = copy_scripts(tdir, testdirs[1]) for fn in files: - run_scripts(fn) + run_tutorial_scripts(fn, testdirs[1]) diff --git a/autotest/t054_test_mfnwt.py b/autotest/t054_test_mfnwt.py index 716e7fcd..0f1988f3 100644 --- a/autotest/t054_test_mfnwt.py +++ b/autotest/t054_test_mfnwt.py @@ -4,6 +4,7 @@ These are the examples that are distributed with MODFLOW-USG. """ import os +import sys import flopy import pymake @@ -37,6 +38,10 @@ v = flopy.which(mfnwt_exe) run = True if v is None: run = False +# fix for intermittent CI failure on windows +else: + if sys.platform.lower() in ("win32", "darwin"): + run = False # diff --git a/examples/Tutorials/Tutorial01/tutorial01.py b/examples/Tutorials/Tutorial01/tutorial01.py deleted file mode 100644 index f08ab8b6..00000000 --- a/examples/Tutorials/Tutorial01/tutorial01.py +++ /dev/null @@ -1,81 +0,0 @@ -import numpy as np -import flopy - -# Assign name and create modflow model object -modelname = "tutorial1" -mf = flopy.modflow.Modflow(modelname, exe_name="mf2005") - -# Model domain and grid definition -Lx = 1000.0 -Ly = 1000.0 -ztop = 0.0 -zbot = -50.0 -nlay = 1 -nrow = 10 -ncol = 10 -delr = Lx / ncol -delc = Ly / nrow -delv = (ztop - zbot) / nlay -botm = np.linspace(ztop, zbot, nlay + 1) - -# Create the discretization object -dis = flopy.modflow.ModflowDis( - mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm[1:] -) - -# Variables for the BAS package -ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) -ibound[:, :, 0] = -1 -ibound[:, :, -1] = -1 -strt = np.ones((nlay, nrow, ncol), dtype=np.float32) -strt[:, :, 0] = 10.0 -strt[:, :, -1] = 0.0 -bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) - -# Add LPF package to the MODFLOW model -lpf = flopy.modflow.ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53) - -# Add OC package to the MODFLOW model -spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} -oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True) - -# Add PCG package to the MODFLOW model -pcg = flopy.modflow.ModflowPcg(mf) - -# Write the MODFLOW model input files -mf.write_input() - -# Run the MODFLOW model -success, buff = mf.run_model() - -# Post process the results -import matplotlib.pyplot as plt -import flopy.utils.binaryfile as bf - -plt.subplot(1, 1, 1, aspect="equal") -hds = bf.HeadFile(modelname + ".hds") -head = hds.get_data(totim=1.0) -levels = np.arange(1, 10, 1) -extent = (delr / 2.0, Lx - delr / 2.0, Ly - delc / 2.0, delc / 2.0) -plt.contour(head[0, :, :], levels=levels, extent=extent) -plt.savefig("tutorial1a.png") - -fig = plt.figure(figsize=(10, 10)) -ax = fig.add_subplot(1, 1, 1, aspect="equal") - -hds = bf.HeadFile(modelname + ".hds") -times = hds.get_times() -head = hds.get_data(totim=times[-1]) -levels = np.linspace(0, 10, 11) - -cbb = bf.CellBudgetFile(modelname + ".cbc") -kstpkper_list = cbb.get_kstpkper() -frf = cbb.get_data(text="FLOW RIGHT FACE", totim=times[-1])[0] -fff = cbb.get_data(text="FLOW FRONT FACE", totim=times[-1])[0] - -modelmap = flopy.plot.ModelMap(model=mf, layer=0) -qm = modelmap.plot_ibound() -lc = modelmap.plot_grid() -cs = modelmap.contour_array(head, levels=levels) -quiver = modelmap.plot_discharge(frf, fff, head=head) -plt.savefig("tutorial1b.png") diff --git a/examples/Tutorials/tutorial01_mf.py b/examples/Tutorials/tutorial01_mf.py new file mode 100644 index 00000000..d36eb6e3 --- /dev/null +++ b/examples/Tutorials/tutorial01_mf.py @@ -0,0 +1,187 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.5.1 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# # MODFLOW Tutorial 1: Confined Steady-State Flow Model +# +# This tutorial demonstrates use of FloPy to develop a simple MODFLOW-2005 +# model. + +# ## Getting Started +# +# If FloPy has been properly installed, then it can be imported as follows: + +import numpy as np +import flopy + +# Now that we can import flopy, we begin creating our simple MODFLOW model. +# numpy is imported to create arrays of model data. + +# ## Creating the MODFLOW Model +# +# One of the nice things about creating models in python is that it is very +# easy to change one or two things and completely change the grid resolution +# for your model. So in this example, we will design our python script so +# that the number of layers, columns, and rows can be easily changed. +# +# We can create a very simple MODFLOW model that has a basic package (BAS), +# discretization input file (DIS), layer-property flow (LPF) package, output +# control (OC), and preconditioned conjugate gradient (PCG) solver. Each one +# of these has its own input file, which will be created automatically by +# flopy, provided that we pass flopy the correct information. + +# ### Discretization +# +# We start by creating our flopy model object. + +modelname = "tutorial1_mf" +mf = flopy.modflow.Modflow(modelname, exe_name="mf2005") + +# Next, let's proceed by defining our model domain and creating a MODFLOW grid +# to span the domain. + +Lx = 1000.0 +Ly = 1000.0 +ztop = 0.0 +zbot = -50.0 +nlay = 1 +nrow = 10 +ncol = 10 +delr = Lx / ncol +delc = Ly / nrow +delv = (ztop - zbot) / nlay +botm = np.linspace(ztop, zbot, nlay + 1) + +# With this information, we can now create the flopy discretization object by +# entering the following: + +dis = flopy.modflow.ModflowDis( + mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm[1:] +) + +# ### Basic Package +# +# Next we can create a flopy object that represents the MODFLOW Basic Package. +# For this simple model, we will assign constant head values of 10. and 0. to +# the first and last model columns (in all layers), respectively. The python +# code for doing this is: + +ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) +ibound[:, :, 0] = -1 +ibound[:, :, -1] = -1 +strt = np.ones((nlay, nrow, ncol), dtype=np.float32) +strt[:, :, 0] = 10.0 +strt[:, :, -1] = 0.0 +bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) + +# ### Layer-Property Flow Package + +# Constant values of 10. are assigned for the horizontal and vertical +# hydraulic conductivity: + +lpf = flopy.modflow.ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53) + +# Because we did not specify a value for laytyp, Flopy will use the default +# value of 0, which means that this model will be confined. + +# ### Output Control +# +# Here we can use the default OC settings by specifying the following: + +spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} +oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True) + +# The stress period dictionary is used to set what output is saved for the +# corresponding stress period and time step. In this case, the tuple `(0, 0)` +# means that stress period 1 and time step 1 for MODFLOW will have output +# saved. Head and budgets will be printed and head and budget information +# will be saved. + +# ### Preconditioned Conjugate Gradient Package +# +# The default settings used by flopy will be used by specifying the following +# commands: + +pcg = flopy.modflow.ModflowPcg(mf) + +# ### Writing the MODFLOW Data Files +# +# The MODFLOW input data files are written by simply issuing the following: + +mf.write_input() + +# ## Running the Model +# +# Flopy can also be used to run the model. The model object (`mf` in this +# example) has an attached method that will run the model. For this to work, +# the MODFLOW program must be located somewhere within the system path, or +# within the working directory. In this example, we have specified that the +# name of the executable program is 'mf2005'. Issue the following to run +# the model: + +success, buff = mf.run_model() +if not success: + raise Exception("MODFLOW did not terminate normally.") + +# Here we have used run_model, and we could also have specified values for the +# optional keywords silent, pause, and report. + +# ## Post-Processing the Results +# +# Now that we have successfully built and run our MODFLOW model, we can look at +# the results. MODFLOW writes the simulated heads to a binary data output file. +# We cannot look at these heads with a text editor, but flopy has a binary +# utility that can be used to read the heads. The following statements will +# read the binary head file and create a plot of simulated heads for layer 1: + +import matplotlib.pyplot as plt +import flopy.utils.binaryfile as bf + +# Extract the heads + +hds = bf.HeadFile(modelname + ".hds") +head = hds.get_data(totim=1.0) + +# Contour the heads + +extent = (delr / 2.0, Lx - delr / 2.0, Ly - delc / 2.0, delc / 2.0) +fig = plt.figure(figsize=(6, 6)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") +ax.contour(head[0, :, :], levels=np.arange(1, 10, 1), extent=extent) + +# Flopy also has some pre-canned plotting capabilities can can be accessed +# using the `PlotMapView()` class. The following code shows how to use the +# plotmapview class to plot boundary conditions (`IBOUND`), plot the grid, +# plot head contours, and plot vectors: + +# Extract the heads +hds = bf.HeadFile(modelname + ".hds") +times = hds.get_times() +head = hds.get_data(totim=times[-1]) + +# Extract the cell-by-cell flows + +cbb = bf.CellBudgetFile(modelname + ".cbc") +kstpkper_list = cbb.get_kstpkper() +frf = cbb.get_data(text="FLOW RIGHT FACE", totim=times[-1])[0] +fff = cbb.get_data(text="FLOW FRONT FACE", totim=times[-1])[0] + +# Create the figure + +fig = plt.figure(figsize=(6, 6)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") +modelmap = flopy.plot.PlotMapView(model=mf, layer=0, ax=ax) +qm = modelmap.plot_ibound() +lc = modelmap.plot_grid() +cs = modelmap.contour_array(head, levels=np.linspace(0, 10, 11)) +quiver = modelmap.plot_discharge(frf, fff, head=head) diff --git a/examples/Tutorials/tutorial01_mf6.py b/examples/Tutorials/tutorial01_mf6.py new file mode 100644 index 00000000..8bd541d7 --- /dev/null +++ b/examples/Tutorials/tutorial01_mf6.py @@ -0,0 +1,255 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.5.1 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# # MODFLOW 6 Tutorial 1: Unconfined Steady-State Flow Model +# +# This tutorial demonstrates use of FloPy to develop a simple MODFLOW 6 +# model. + +# ## Getting Started + +import os +import numpy as np +import matplotlib.pyplot as plt +import flopy + +# We are creating a square model with a specified head equal to `h1` along +# all boundaries. The head at the cell in the center in the top layer is +# fixed to `h2`. First, set the name of the model and the parameters of the +# model: the number of layers `Nlay`, the number of rows and columns `N`, +# lengths of the sides of the model `L`, aquifer thickness `H`, +# hydraulic conductivity `k` + +name = "tutorial01_mf6" +h1 = 100 +h2 = 90 +Nlay = 10 +N = 101 +L = 400.0 +H = 50.0 +k = 1.0 + +# ### Create the Flopy Model Objects +# +# One big difference between MODFLOW 6 and previous MODFLOW versions is +# that MODFLOW 6 is based on the concept of a simulation. +# A simulation consists of the following: +# +# * Temporal discretization (`TDIS`) +# * One or more models +# * Zero or more exchanges (instructions for how models are coupled) +# * Solutions +# +# For this simple example, the simulation consists of the temporal +# discretization (`TDIS`) package, a groundwater flow (`GWF`) model, and +# an iterative model solution (`IMS`), which controls how the GWF model is +# solved. + +# ### Create the Flopy simulation object + +sim = flopy.mf6.MFSimulation( + sim_name=name, exe_name="mf6", version="mf6", sim_ws="." +) + +# ### Create the Flopy `TDIS` object + +tdis = flopy.mf6.ModflowTdis( + sim, pname="tdis", time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)] +) + +# ### Create the Flopy `IMS` Package object + +ims = flopy.mf6.ModflowIms(sim, pname="ims", complexity="SIMPLE") + +# Create the Flopy groundwater flow (gwf) model object + +model_nam_file = "{}.nam".format(name) +gwf = flopy.mf6.ModflowGwf(sim, modelname=name, model_nam_file=model_nam_file) + +# Now that the overall simulation is set up, we can focus on building the +# groundwater flow model. The groundwater flow model will be built by +# adding packages to it that describe the model characteristics. + +# ### Create the discretization (`DIS`) Package +# +# Define the discretization of the model. All layers are given equal thickness. +# The `bot` array is build from `H` and the `Nlay` values to indicate top and +# bottom of each layer, and `delrow` and `delcol` are computed from model +# size `L` and number of cells `N`. Once these are all computed, the +# Discretization file is built. + +bot = np.linspace(-H / Nlay, -H, Nlay) +delrow = delcol = L / (N - 1) +dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=Nlay, + nrow=N, + ncol=N, + delr=delrow, + delc=delcol, + top=0.0, + botm=bot, +) + +# ### Create the initial conditions (`IC`) Package + +start = h1 * np.ones((Nlay, N, N)) +ic = flopy.mf6.ModflowGwfic(gwf, pname="ic", strt=start) + +# ### Create the node property flow (`NPF`) Package + +npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=k, save_flows=True) + +# ### Create the constant head (`CHD`) Package +# +# List information is created a bit differently for MODFLOW 6 than for other +# MODFLOW versions. The cellid (layer, row, column, for a regular grid) +# can be entered as a tuple as the first entry. Remember that these must be +# zero-based indices! + +chd_rec = [] +chd_rec.append(((0, int(N / 4), int(N / 4)), h2)) +for layer in range(0, Nlay): + for row_col in range(0, N): + chd_rec.append(((layer, row_col, 0), h1)) + chd_rec.append(((layer, row_col, N - 1), h1)) + if row_col != 0 and row_col != N - 1: + chd_rec.append(((layer, 0, row_col), h1)) + chd_rec.append(((layer, N - 1, row_col), h1)) +chd = flopy.mf6.ModflowGwfchd( + gwf, + maxbound=len(chd_rec), + stress_period_data=chd_rec, + save_flows=True, +) + +# The `CHD` Package stored the constant heads in a structured array, +# also called a `numpy.recarray`. We can get a pointer to the recarray +# for the first stress period (iper = 0) as follows. + +iper = 0 +ra = chd.stress_period_data.get_data(key=iper) +ra + +# Create the output control (`OC`) Package +headfile = "{}.hds".format(name) +head_filerecord = [headfile] +budgetfile = "{}.cbb".format(name) +budget_filerecord = [budgetfile] +saverecord = [("HEAD", "ALL"), ("BUDGET", "ALL")] +printrecord = [("HEAD", "LAST")] +oc = flopy.mf6.ModflowGwfoc( + gwf, + saverecord=saverecord, + head_filerecord=head_filerecord, + budget_filerecord=budget_filerecord, + printrecord=printrecord, +) + + +# ## Create the MODFLOW 6 Input Files and Run the Model +# +# Once all the flopy objects are created, it is very easy to create +# all of the input files and run the model. + + +# ### Write the datasets + +sim.write_simulation() + +# ### Run the Simulation +# +# We can also run the simulation from python, but only if the MODFLOW 6 +# executable is available. The executable can be made available by putting +# the executable in a folder that is listed in the system path variable. +# Another option is to just put a copy of the executable in the simulation +# folder, though this should generally be avoided. A final option is to +# provide a full path to the executable when the simulation is constructed. +# This would be done by specifying exe_name with the full path. + +success, buff = sim.run_simulation() +if not success: + raise Exception("MODFLOW 6 did not terminate normally.") + +# ## Post-Process Head Results +# +# First, a link to the heads file is created with `HeadFile`. The link can then be accessed with the `get_data` function, by specifying, in this case, the step number and period number for which we want to retrieve data. A three-dimensional array is returned of size `nlay, nrow, ncol`. Matplotlib contouring functions are used to make contours of the layers or a cross-section. + +# Read the binary head file and plot the results. We can use the Flopy +# `HeadFile()` class because the format of the headfile for MODFLOW 6 is +# the same as for previous MODFLOW versions. + +# ### Plot a Map of Layer 1 + +hds = flopy.utils.binaryfile.HeadFile(headfile) +h = hds.get_data(kstpkper=(0, 0)) +x = y = np.linspace(0, L, N) +y = y[::-1] +fig = plt.figure(figsize=(6, 6)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") +c = ax.contour(x, y, h[0], np.arange(90, 100.1, 0.2), colors="black") +plt.clabel(c, fmt="%2.1f") + + +# ### Plot a Map of Layer 10 + +x = y = np.linspace(0, L, N) +y = y[::-1] +fig = plt.figure(figsize=(6, 6)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") +c = ax.contour(x, y, h[-1], np.arange(90, 100.1, 0.2), colors="black") +plt.clabel(c, fmt="%1.1f") + +# ### Plot a Cross-section along row 51 + +z = np.linspace(-H / Nlay / 2, -H + H / Nlay / 2, Nlay) +fig = plt.figure(figsize=(5, 2.5)) +ax = fig.add_subplot(1, 1, 1, aspect="auto") +c = ax.contour(x, z, h[:, 50, :], np.arange(90, 100.1, 0.2), colors="black") +plt.clabel(c, fmt="%1.1f") + +# ### We can also use the Flopy `PlotMapView()` capabilities for MODFLOW 6 +# +# Before we start we will create a MODFLOW-2005 ibound array to use to plot +# the locations of the constant heads. + +ibd = np.ones((Nlay, N, N), dtype=np.int) +for k, i, j in ra["cellid"]: + ibd[k, i, j] = -1 + +# ### Plot a Map of Layers 1 and 10 + +fig, axes = plt.subplots(2, 1, figsize=(6, 12), constrained_layout=True) +# first subplot +ax = axes[0] +ax.set_title("Model Layer 1") +modelmap = flopy.plot.PlotMapView(model=gwf, ax=ax) +quadmesh = modelmap.plot_ibound(ibound=ibd) +linecollection = modelmap.plot_grid(lw=0.5, color="0.5") +contours = modelmap.contour_array( + h[0], levels=np.arange(90, 100.1, 0.2), colors="black" +) +ax.clabel(contours, fmt="%2.1f") +# second subplot +ax = axes[1] +ax.set_title("Model Layer {}".format(Nlay)) +modelmap = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=Nlay - 1) +quadmesh = modelmap.plot_ibound(ibound=ibd) +linecollection = modelmap.plot_grid(lw=0.5, color="0.5") +pa = modelmap.plot_array(h[0]) +contours = modelmap.contour_array( + h[0], levels=np.arange(90, 100.1, 0.2), colors="black" +) +cb = plt.colorbar(pa, shrink=0.5, ax=ax) +ax.clabel(contours, fmt="%2.1f") diff --git a/examples/Tutorials/tutorial01_seawat.py b/examples/Tutorials/tutorial01_seawat.py new file mode 100644 index 00000000..e33698bd --- /dev/null +++ b/examples/Tutorials/tutorial01_seawat.py @@ -0,0 +1,196 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.5.1 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# # SEAWAT Tutorial 1: Henry Saltwater Intrusion Problem +# +# In this tutorial, we will use Flopy to create, run, and post process the +# Henry saltwater intrusion problem using SEAWAT Version 4. + +# ## Getting Started + +import os +import numpy as np +import flopy + +# ### Input variables for the Henry Problem + +Lx = 2.0 +Lz = 1.0 +nlay = 50 +nrow = 1 +ncol = 100 +delr = Lx / ncol +delc = 1.0 +delv = Lz / nlay +henry_top = 1.0 +henry_botm = np.linspace(henry_top - delv, 0.0, nlay) +qinflow = 5.702 # m3/day +dmcoef = 0.57024 # m2/day Could also try 1.62925 as another case of the Henry problem +hk = 864.0 # m/day + +# ### Create the basic MODFLOW model structure + +modelname = "henry" +swt = flopy.seawat.Seawat(modelname, exe_name="swtv4") +print(swt.namefile) + +# save cell fluxes to unit 53 + +ipakcb = 53 + +# Add DIS package to the MODFLOW model + +dis = flopy.modflow.ModflowDis( + swt, + nlay, + nrow, + ncol, + nper=1, + delr=delr, + delc=delc, + laycbd=0, + top=henry_top, + botm=henry_botm, + perlen=1.5, + nstp=15, +) + +# Variables for the BAS package +ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) +ibound[:, :, -1] = -1 + +# Add BAS package to the MODFLOW model + +bas = flopy.modflow.ModflowBas(swt, ibound, 0) + +# Add LPF package to the MODFLOW model + +lpf = flopy.modflow.ModflowLpf(swt, hk=hk, vka=hk, ipakcb=ipakcb) + +# Add PCG Package to the MODFLOW model + +pcg = flopy.modflow.ModflowPcg(swt, hclose=1.0e-8) + +# Add OC package to the MODFLOW model + +oc = flopy.modflow.ModflowOc( + swt, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, +) + +# Create WEL and SSM data + +itype = flopy.mt3d.Mt3dSsm.itype_dict() +wel_data = {} +ssm_data = {} +wel_sp1 = [] +ssm_sp1 = [] +for k in range(nlay): + wel_sp1.append([k, 0, 0, qinflow / nlay]) + ssm_sp1.append([k, 0, 0, 0.0, itype["WEL"]]) + ssm_sp1.append([k, 0, ncol - 1, 35.0, itype["BAS6"]]) +wel_data[0] = wel_sp1 +ssm_data[0] = ssm_sp1 +wel = flopy.modflow.ModflowWel(swt, stress_period_data=wel_data, ipakcb=ipakcb) + + +# ### Create the basic MT3DMS model structure + +btn = flopy.mt3d.Mt3dBtn( + swt, + nprs=-5, + prsity=0.35, + sconc=35.0, + ifmtcn=0, + chkmas=False, + nprobs=10, + nprmas=10, + dt0=0.001, +) +adv = flopy.mt3d.Mt3dAdv(swt, mixelm=0) +dsp = flopy.mt3d.Mt3dDsp(swt, al=0.0, trpt=1.0, trpv=1.0, dmcoef=dmcoef) +gcg = flopy.mt3d.Mt3dGcg(swt, iter1=500, mxiter=1, isolve=1, cclose=1e-7) +ssm = flopy.mt3d.Mt3dSsm(swt, stress_period_data=ssm_data) + +# ### Create the SEAWAT model structure + +vdf = flopy.seawat.SeawatVdf( + swt, + iwtable=0, + densemin=0, + densemax=0, + denseref=1000.0, + denseslp=0.7143, + firstdt=1e-3, +) + +# ## Write the input files + +swt.write_input() + +# ## Run the model + +success, buff = swt.run_model(silent=True, report=True) +if not success: + raise Exception("SEAWAT did not terminate normally.") + +# ## Post-process the results + +import numpy as np +import flopy.utils.binaryfile as bf + +# ### Load the concentration data + +ucnobj = bf.UcnFile("MT3D001.UCN", model=swt) +times = ucnobj.get_times() +concentration = ucnobj.get_data(totim=times[-1]) + +# ### Load the cell-by-cell flow data + +cbbobj = bf.CellBudgetFile("henry.cbc") +times = cbbobj.get_times() +qx = cbbobj.get_data(text="flow right face", totim=times[-1])[0] +qy = np.zeros((nlay, nrow, ncol), dtype=np.float) +qz = cbbobj.get_data(text="flow lower face", totim=times[-1])[0] + +# ### Create a plot with concentrations and flow vectors + +import matplotlib.pyplot as plt + +fig = plt.figure(figsize=(12,9)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") +pmv = flopy.plot.PlotCrossSection(model=swt, ax=ax, line={"row": 0}) +arr = pmv.plot_array(concentration) +pmv.plot_vector(qx, qy, -qz, color="white", kstep=3, hstep=3) +plt.colorbar(arr, shrink=0.5, ax=ax) +ax.set_title("Simulated Concentrations"); + +# ### Load the head data + +headobj = bf.HeadFile("henry.hds") +times = headobj.get_times() +head = headobj.get_data(totim=times[-1]) + + +# ### Create a plot with heads + +fig = plt.figure(figsize=(12, 9)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") +pmv = flopy.plot.PlotCrossSection(model=swt, ax=ax, line={"row": 0}) +arr = pmv.plot_array(head) +contours = pmv.contour_array(head, colors="white") +ax.clabel(contours, fmt="%2.2f") +plt.colorbar(arr, shrink=0.5, ax=ax) +ax.set_title("Simulated Heads"); diff --git a/examples/Tutorials/Tutorial02/tutorial02.py b/examples/Tutorials/tutorial02_mf.py similarity index 50% rename from examples/Tutorials/Tutorial02/tutorial02.py rename to examples/Tutorials/tutorial02_mf.py index af0c8378..ab2cfd44 100644 --- a/examples/Tutorials/Tutorial02/tutorial02.py +++ b/examples/Tutorials/tutorial02_mf.py @@ -1,7 +1,48 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.5.1 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# # MODFLOW Tutorial 2: Unconfined Transient Flow Model +# +# In this example, we will convert the tutorial 1 model into an unconfined, +# transient flow model with time varying boundaries. Instead of using constant +# heads for the left and right boundaries (by setting ibound to -1), we will use +# general head boundaries. We will have the model consider the following +# conditions: +# +# * Initial conditions -- head is 10.0 everywhere +# * Period 1 (1 day) -- steady state with left and right GHB stage = 10. +# * Period 2 (100 days) -- left GHB with stage = 10., right GHB with stage set +# to 0. +# * Period 3 (100 days) -- pumping well at model center with rate = -500., left +# and right GHB = 10., and 0. +# +# We will start with selected model commands from the previous tutorial. +# + +# ## Getting Started +# +# As shown in the previous MODFLOW tutorial, import flopy. + import numpy as np import flopy -# Model domain and grid definition +# ## Creating the MODFLOW Model +# +# ### Define the Model Extent, Grid Resolution, and Characteristics +# +# Assign the model information + Lx = 1000.0 Ly = 1000.0 ztop = 10.0 @@ -20,18 +61,28 @@ ss = 1.0e-4 laytyp = 1 # Variables for the BAS package -# Note that changes from the previous tutorial! +# Note that changes from the MODFLOW tutorial 1 + ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) strt = 10.0 * np.ones((nlay, nrow, ncol), dtype=np.float32) -# Time step parameters +# ### Define the Stress Periods +# +# To create a model with multiple stress periods, we need to define nper, +# perlen, nstp, and steady. This is done in the following block in a manner +# that allows us to pass these variable directly to the discretization object: + nper = 3 perlen = [1, 100, 100] nstp = [1, 100, 100] steady = [True, False, False] -# Flopy objects -modelname = "tutorial2" +# ### Create Time-Invariant Flopy Objects +# +# With this information, we can now create the static flopy objects that do +# not change with time: + +modelname = "tutorial2_mf" mf = flopy.modflow.Modflow(modelname, exe_name="mf2005") dis = flopy.modflow.ModflowDis( mf, @@ -53,7 +104,18 @@ lpf = flopy.modflow.ModflowLpf( ) pcg = flopy.modflow.ModflowPcg(mf) -# Make list for stress period 1 +# ### Transient General-Head Boundary Package +# +# At this point, our model is ready to add our transient boundary packages. +# First, we will create the GHB object, which is of the following type: +# `flopy.modflow.ModflowGhb()`. +# +# The key to creating Flopy transient boundary packages is recognizing that +# the boundary data is stored in a dictionary with key values equal to the +# zero-based stress period number and values equal to the boundary conditions +# for that stress period. For a GHB the values can be a two-dimensional nested +# list of `[layer, row, column, stage, conductance]`. + stageleft = 10.0 stageright = 10.0 bound_sp1 = [] @@ -85,6 +147,11 @@ stress_period_data = {0: bound_sp1, 1: bound_sp2} # Create the flopy ghb object ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=stress_period_data) +# ### Transient Well Package +# +# Now we can create the well package object, which is of the type, +# `flopy.modflow.ModflowWel()`. + # Create the well package # Remember to use zero-based layer, row, column indices! pumping_rate = -500.0 @@ -94,7 +161,11 @@ wel_sp3 = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]] stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3} wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data) -# Output control +# ### Output Control +# +# Here we create the output control package object, which is of the +# type `flopy.modflow.ModflowOc()`. + stress_period_data = {} for kper in range(nper): for kstp in range(nstp[kper]): @@ -109,14 +180,34 @@ oc = flopy.modflow.ModflowOc( mf, stress_period_data=stress_period_data, compact=True ) +# ## Running the Model +# +# Run the model with the run_model method, which returns a success flag and +# the stream of output. With run_model, we have some finer control, that +# allows us to suppress the output. + # Write the model input files mf.write_input() # Run the model -success, mfoutput = mf.run_model(silent=False, pause=False) +success, mfoutput = mf.run_model(silent=True, pause=False) if not success: raise Exception("MODFLOW did not terminate normally.") +# ## Post-Processing the Results +# +# Once again, we can read heads from the MODFLOW binary output file, using +# the `flopy.utils.binaryfile()` module. Included with the HeadFile object +# are several methods that we will use here: + +# * `get_times()` will return a list of times contained in the binary head file +# * `get_data()` will return a three-dimensional head array for the specified +# time +# * `get_ts()` will return a time series array `[ntimes, headval]` for the +# specified cell +# +# Using these methods, we can create head plots and hydrographs from the +# model results. # Imports import matplotlib.pyplot as plt @@ -133,11 +224,13 @@ extent = (delr / 2.0, Lx - delr / 2.0, delc / 2.0, Ly - delc / 2.0) print("Levels: ", levels) print("Extent: ", extent) -# Well point -wpt = ((float(ncol / 2) - 0.5) * delr, (float(nrow / 2 - 1) + 0.5) * delc) +# Well point for plotting wpt = (450.0, 550.0) +# Create a figure with maps for three times + # Make the plots +fig = plt.figure(figsize=(5, 15)) mytimes = [1.0, 101.0, 201.0] for iplot, time in enumerate(mytimes): print("*****Processing time: ", time) @@ -152,23 +245,23 @@ for iplot, time in enumerate(mytimes): frf = cbb.get_data(text="FLOW RIGHT FACE", totim=time)[0] fff = cbb.get_data(text="FLOW FRONT FACE", totim=time)[0] - # Create the plot - # plt.subplot(1, len(mytimes), iplot + 1, aspect='equal') - plt.subplot(1, 1, 1, aspect="equal") - plt.title("stress period " + str(iplot + 1)) + # Create a map for this time + ax = fig.add_subplot(len(mytimes), 1, iplot + 1, aspect="equal") + ax.set_title("stress period " + str(iplot + 1)) - modelmap = flopy.plot.ModelMap(model=mf, layer=0) - qm = modelmap.plot_ibound() - lc = modelmap.plot_grid() - qm = modelmap.plot_bc("GHB", alpha=0.5) - cs = modelmap.contour_array(head, levels=levels) - plt.clabel(cs, inline=1, fontsize=10, fmt="%1.1f") - quiver = modelmap.plot_discharge(frf, fff, head=head) + pmv = flopy.plot.PlotMapView(model=mf, layer=0, ax=ax) + qm = pmv.plot_ibound() + lc = pmv.plot_grid() + qm = pmv.plot_bc("GHB", alpha=0.5) + if head.min() != head.max(): + cs = pmv.contour_array(head, levels=levels) + plt.clabel(cs, inline=1, fontsize=10, fmt="%1.1f") + quiver = pmv.plot_vector(frf, fff) mfc = "None" if (iplot + 1) == len(mytimes): mfc = "black" - plt.plot( + ax.plot( wpt[0], wpt[1], lw=0, @@ -179,17 +272,17 @@ for iplot, time in enumerate(mytimes): markerfacecolor=mfc, zorder=9, ) - plt.text(wpt[0] + 25, wpt[1] - 25, "well", size=12, zorder=12) - plt.savefig("tutorial2-{}.png".format(iplot)) + ax.text(wpt[0] + 25, wpt[1] - 25, "well", size=12, zorder=12) +# Create a hydrograph # Plot the head versus time idx = (0, int(nrow / 2) - 1, int(ncol / 2) - 1) ts = headobj.get_ts(idx) -plt.subplot(1, 1, 1) +fig = plt.figure(figsize=(6, 6)) +ax = fig.add_subplot(1, 1, 1) ttl = "Head at cell ({0},{1},{2})".format(idx[0] + 1, idx[1] + 1, idx[2] + 1) -plt.title(ttl) -plt.xlabel("time") -plt.ylabel("head") -plt.plot(ts[:, 0], ts[:, 1], "bo-") -plt.savefig("tutorial2-ts.png") +ax.set_title(ttl) +ax.set_xlabel("time") +ax.set_ylabel("head") +ax.plot(ts[:, 0], ts[:, 1], "bo-")