From 510b9f169f35e60f6850ea2a66df0aed277a441d Mon Sep 17 00:00:00 2001 From: Chris Nicol <37314969+cnicol-gwlogic@users.noreply.github.com> Date: Mon, 10 Aug 2020 22:45:38 +1000 Subject: [PATCH] refactor(gridgen): fixes and additions related to vertical pass through cells (#953) * Gridgen fixes and additions - Added Gridgen's "vertical pass through" option to utils/gridgen.export(). - Fixed utils/gridgen.set_surface_interpolation() file paths for ASCIIGRID option. - Fixed grd variable definitions for ASCIIGRID surface interpolation option in utils/gridgen._builder_block(). These now use gridgen._asciigrid_dict where appropriate (default ascii grid names were causing errors). - Added docstring to utils/gridgen.get_disu(). - Minor edit to modflow/mfdisu.py docstring re: util3d / unstructured models. * docs(readme): autotest build_exes.py changed to get_exes.py * test(gridgen vertical pass through) * test(gridgen vertical pass through) * style(gridgen.py) * test(gridgen vertical pass through) --- CONTRIBUTING.md | 2 +- autotest/t061_test_gridgen.py | 49 ++++++++++++++++++++++++++++ flopy/modflow/mfdisu.py | 5 ++- flopy/utils/gridgen.py | 61 ++++++++++++++++++++++++++++++++--- 4 files changed, 108 insertions(+), 9 deletions(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 0bf2bf77..f0ce6a1d 100755 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -66,7 +66,7 @@ Before you submit your Pull Request (PR) consider the following guidelines: ```shell cd autotest - nosetests -v build_exes.py + nosetests -v get_exes.py nosetests -v ``` Note: the FloPy test suite requires the [nosetests](https://pypi.org/project/nose/) and [pymake](https://github.com/modflowpy/pymake) python packages. All the FloPy dependencies must also be installed for the tests to pass. diff --git a/autotest/t061_test_gridgen.py b/autotest/t061_test_gridgen.py index b2b1c88f..2d067fb7 100644 --- a/autotest/t061_test_gridgen.py +++ b/autotest/t061_test_gridgen.py @@ -45,9 +45,17 @@ def test_gridgen(): delr=delr, delc=delc, top=top, botm=botm) + ms_u = flopy.modflow.Modflow(modelname = 'mymfusgmodel', model_ws = cpth, + version = 'mfusg') + dis_usg = flopy.modflow.ModflowDis(ms_u, nlay=nlay, nrow=nrow, ncol=ncol, + delr=delr, + delc=delc, top=top, botm=botm) + gridgen_ws = cpth g = Gridgen(dis5, model_ws=gridgen_ws, exe_name=exe_name) g6 = Gridgen(dis6, model_ws=gridgen_ws, exe_name=exe_name) + gu = Gridgen(dis_usg, model_ws=gridgen_ws, exe_name=exe_name, + vertical_pass_through=True) rf0shp = os.path.join(gridgen_ws, 'rf0') xmin = 7 * delr @@ -58,6 +66,7 @@ def test_gridgen(): (xmin, ymin)]]] g.add_refinement_features(rfpoly, 'polygon', 1, range(nlay)) g6.add_refinement_features(rfpoly, 'polygon', 1, range(nlay)) + gu.add_refinement_features(rfpoly, 'polygon', 1, range(nlay)) rf1shp = os.path.join(gridgen_ws, 'rf1') xmin = 8 * delr @@ -68,6 +77,7 @@ def test_gridgen(): (xmin, ymin)]]] g.add_refinement_features(rfpoly, 'polygon', 2, range(nlay)) g6.add_refinement_features(rfpoly, 'polygon', 2, range(nlay)) + gu.add_refinement_features(rfpoly, 'polygon', 2, range(nlay)) rf2shp = os.path.join(gridgen_ws, 'rf2') xmin = 9 * delr @@ -78,6 +88,19 @@ def test_gridgen(): (xmin, ymin)]]] g.add_refinement_features(rfpoly, 'polygon', 3, range(nlay)) g6.add_refinement_features(rfpoly, 'polygon', 3, range(nlay)) + gu.add_refinement_features(rfpoly, 'polygon', 3, range(nlay)) + + # inactivate parts of mfusg layer 2 to test vertical-pass-through option + xmin = 0 * delr + xmax = 18 * delr + ymin = 0 * delc + ymax = 18 * delc + adpoly2 = [[[(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax), + (xmin, ymin)]]] + gu.add_active_domain(adpoly2, layers = [1]) + adpoly1_3 = [[[(0., 0.), (Lx, 0.), (Lx, Ly), (0., Ly), + (0., 0.)]]] + gu.add_active_domain(adpoly1_3, layers = [0,2]) # if gridgen executable is available then do the main part of the test if run: @@ -126,6 +149,32 @@ def test_gridgen(): assert os.path.isfile(fname), \ 'MF6 disv file not created: {}'.format(fname) + # test mfusg with vertical pass-through (True above at instantiation) + gu.build() + disu_vp = gu.get_disu(ms_u) + # -check that node 1 (layer 1) is connected to layer 3 but not layer 2: + ja0 = disu_vp.ja[: disu_vp.iac[0]] + msg = ("MFUSG node 1 (layer 1) is not connected to layer 3 but should " + "be (with vertical pass through activated).") + assert max(ja0) > sum(disu_vp.nodelay[:2]), msg + # -check that node 1 (layer 1) is not connected to any layer 2 nodes + msg = ("MFUSG node 1 (layer 1) is connected to layer 2 but should not " + "be (with vertical pass through activated).") + assert len(ja0[(ja0 > disu_vp.nodelay[0]) & \ + (ja0 <= sum(disu_vp.nodelay[:2]))] + ) == 0, msg + #ms_u.disu.write_file() + + # test mfusg without vertical pass-through + gu.vertical_pass_through = False + gu.build() + disu_vp = gu.get_disu(ms_u) + # -check that node 1 (layer 1) is connected to layer 1 only: + ja0 = disu_vp.ja[: disu_vp.iac[0]] + msg = ("MFUSG node 1 (layer 1) is connected to layer 2 or 3 but " + "should not be (without vertical pass through activated).") + assert max(ja0) <= disu_vp.nodelay[0], msg + return diff --git a/flopy/modflow/mfdisu.py b/flopy/modflow/mfdisu.py index f9179c7a..917fa70b 100644 --- a/flopy/modflow/mfdisu.py +++ b/flopy/modflow/mfdisu.py @@ -184,9 +184,8 @@ class ModflowDisU(Package): Notes ----- - Does not work yet for multi-layer USG models because top and bot cannot - be u3d instances until u3d is modified to handle multiple u2d instances - of different size. + Now works for multi-layer USG models since u3d was modified to handle + multiple u2d instances of different size. Examples -------- diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py index dab283ec..956d1762 100644 --- a/flopy/utils/gridgen.py +++ b/flopy/utils/gridgen.py @@ -121,6 +121,10 @@ class Gridgen(object): surface_interpolation : str Default gridgen method for interpolating elevations. Valid options include 'replicate' (default) and 'interpolate' + vertical_pass_through : bool + If true, Gridgen's GRID_TO_USGDATA command will connect layers + where intermediate layers are inactive. + (default is False) Notes ----- @@ -135,6 +139,7 @@ class Gridgen(object): model_ws=".", exe_name="gridgen", surface_interpolation="replicate", + vertical_pass_through=False, ): self.dis = dis if isinstance(dis, ModflowGwfdis): @@ -170,6 +175,11 @@ class Gridgen(object): surface_interpolation for k in range(self.nlay + 1) ] + # Set export options + self.vertical_pass_through = "False" + if vertical_pass_through: + self.vertical_pass_through = "True" + # Set up a blank _active_domain list with None for each layer self._addict = {} self._active_domain = [] @@ -202,7 +212,8 @@ class Gridgen(object): Array that is used as an asciigrid. If elev is a string, then it is assumed to be the name of the asciigrid. elev_extent : list-like - list of xmin, xmax, ymin, ymax extents of the elev grid. + List of xmin, xmax, ymin, ymax extents of the elev grid. + Must be specified for ASCIIGRID; optional otherwise. Returns ------- @@ -242,9 +253,10 @@ class Gridgen(object): self._asciigrid_dict[isurf] = nm elif isinstance(elev, str): - if not os.path.isfile(elev): + if not os.path.isfile(os.path.join(self.model_ws, elev)): raise Exception( - "Error. elev is not a valid file: " "{}".format(elev) + "Error. elev is not a valid file: " + "{}".format(os.path.join(self.model_ws, elev)) ) self._asciigrid_dict[isurf] = elev else: @@ -439,11 +451,18 @@ class Gridgen(object): """ Export the quadtree grid to shapefiles, usgdata, and vtk + Parameters + ---------- + verbose : bool + If true, print the results of the gridgen command to the terminal + (default is False) + Returns ------- None """ + # Create the export definition file fname = os.path.join(self.model_ws, "_gridgen_export.dfn") f = open(fname, "w") @@ -635,6 +654,35 @@ class Gridgen(object): itmuni=4, lenuni=2, ): + """ + Create a MODFLOW-USG DISU flopy object. + + Parameters + ---------- + model : Flopy model object + The Flopy model object (of type :class:`flopy.modflow.mf.Modflow`) + to which this package will be added. + nper : int + Number of model stress periods (default is 1). + perlen : float or array of floats (nper) + A single value or array of the stress period lengths + (default is 1). + nstp : int or array of ints (nper) + Number of time steps in each stress period (default is 1). + tsmult : float or array of floats (nper) + Time step multiplier (default is 1.0). + steady : boolean or array of boolean (nper) + True or False indicating whether or not stress period is + steady state (default is True). + itmuni : int + Time units, default is days (4) + lenuni : int + Length units, default is meters (2) + + Returns + ------- + disu : Flopy ModflowDisU object. + """ # nodes, nlay, ivsd, itmuni, lenuni, idsymrd, laycbd fname = os.path.join(self.model_ws, "qtg.nod") @@ -1787,7 +1835,7 @@ class Gridgen(object): for k in range(self.nlay): if self.surface_interpolation[k] == "ASCIIGRID": - grd = "_gridgen.lay{}.asc".format(k) + grd = self._asciigrid_dict[k] else: grd = "basename" s += " TOP LAYER {} = {} {}\n".format( @@ -1796,7 +1844,7 @@ class Gridgen(object): for k in range(self.nlay): if self.surface_interpolation[k + 1] == "ASCIIGRID": - grd = "_gridgen.lay{}.asc".format(k + 1) + grd = self._asciigrid_dict[k + 1] else: grd = "basename" s += " BOTTOM LAYER {} = {} {}\n".format( @@ -1823,6 +1871,9 @@ class Gridgen(object): s += "BEGIN GRID_TO_USGDATA grid_to_usgdata\n" s += " GRID = quadtreegrid\n" s += " USG_DATA_PREFIX = qtg\n" + s += " VERTICAL_PASS_THROUGH = {0}\n".format( + self.vertical_pass_through + ) s += "END GRID_TO_USGDATA\n" s += "\n" s += "BEGIN GRID_TO_VTKFILE grid_to_vtk\n"