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)
develop
Chris Nicol 2020-08-10 22:45:38 +10:00 committed by GitHub
parent 05ce3bc53d
commit 510b9f169f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 108 additions and 9 deletions

View File

@ -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.

View File

@ -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

View File

@ -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
--------

View File

@ -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"