Merge pull request #2305 from modflowpy/v3.8.1

Release 3.8.1

Author: wpbonelli

.docs/md/version_changes.md

@@ -1,4 +1,19 @@
 # Changelog
+### Version 3.8.1
+
+#### New features
+
+* [feat(cell1d)](https://github.com/modflowpy/flopy/commit/4ea71927f251c3675acf1b578bca623d023ccd2c): Add support for 1D vertex grids (#2296). Committed by langevin-usgs on 2024-08-23.
+
+#### Bug fixes
+
+* [fix(ParticleTrackFile.write_shapefile)](https://github.com/modflowpy/flopy/commit/f86881d6354071bf7c675384c2684ea184e281eb): Check for "k" even if "i", "j are not present   (#2294). Committed by Joshua Larsen on 2024-08-17.
+* [fix(modelgrid)](https://github.com/modflowpy/flopy/commit/c42d8787bbcd4131b2f493ebc5a5a384b2e8b861): Add more support for mf6 surface water models  (#2295). Committed by langevin-usgs on 2024-08-22.
+
+#### Refactoring
+
+* [refactor(model_splitter.py)](https://github.com/modflowpy/flopy/commit/d02967db167b98e32cfaa26ef6636475ea2441a8): Update UnstructuredGrid support (#2292). Committed by Joshua Larsen on 2024-08-16.
+
 ### Version 3.8.0
 
 #### New features

.github/workflows/release.yml

@@ -84,7 +84,7 @@ jobs:
 
       - name: Generate changelog
         id: cliff
-        uses: orhun/git-cliff-action@v3
+        uses: orhun/git-cliff-action@v4
         with:
           config: cliff.toml
           args: --verbose --unreleased --tag ${{ steps.version.outputs.version }}
@@ -120,7 +120,7 @@ jobs:
           git config core.sharedRepository true
           git config user.name "github-actions[bot]"
           git config user.email "41898282+github-actions[bot]@users.noreply.github.com"
-          git add -A
+          git add flopy docs .docs CITATION.cff README.md version.txt 
           git commit -m "ci(release): set version to ${{ steps.version.outputs.version }}, update plugins from DFN files, update changelog"
           git push origin "${{ github.ref_name }}"
 

CITATION.cff

@@ -3,8 +3,8 @@ message: If you use this software, please cite both the article from preferred-c
   references, and the software itself.
 type: software
 title: FloPy
-version: 3.8.0
-date-released: '2024-08-08'
+version: 3.8.1
+date-released: '2024-09-05'
 doi: 10.5066/F7BK19FH
 abstract: A Python package to create, run, and post-process MODFLOW-based models.
 repository-artifact: https://pypi.org/project/flopy

README.md

@@ -1,7 +1,7 @@
 
 <img src="https://raw.githubusercontent.com/modflowpy/flopy/master/examples/images/flopy3.png" alt="flopy3" style="width:50;height:20">
 
-### Version 3.8.0
+### Version 3.8.1
 [![flopy continuous integration](https://github.com/modflowpy/flopy/actions/workflows/commit.yml/badge.svg?branch=develop)](https://github.com/modflowpy/flopy/actions/workflows/commit.yml)
 [![Read the Docs](https://github.com/modflowpy/flopy/actions/workflows/rtd.yml/badge.svg?branch=develop)](https://github.com/modflowpy/flopy/actions/workflows/rtd.yml)
 
@@ -150,7 +150,7 @@ How to Cite
 
 ##### ***Software/Code citation for FloPy:***
 
-[Bakker, Mark, Post, Vincent, Hughes, J. D., Langevin, C. D., White, J. T., Leaf, A. T., Paulinski, S. R., Bellino, J. C., Morway, E. D., Toews, M. W., Larsen, J. D., Fienen, M. N., Starn, J. J., Brakenhoff, D. A., and Bonelli, W. P., 2024, FloPy v3.8.0: U.S. Geological Survey Software Release, 08 August 2024, https://doi.org/10.5066/F7BK19FH](https://doi.org/10.5066/F7BK19FH)
+[Bakker, Mark, Post, Vincent, Hughes, J. D., Langevin, C. D., White, J. T., Leaf, A. T., Paulinski, S. R., Bellino, J. C., Morway, E. D., Toews, M. W., Larsen, J. D., Fienen, M. N., Starn, J. J., Brakenhoff, D. A., and Bonelli, W. P., 2024, FloPy v3.8.1: U.S. Geological Survey Software Release, 05 September 2024, https://doi.org/10.5066/F7BK19FH](https://doi.org/10.5066/F7BK19FH)
 
 
 Additional FloPy Related Publications

autotest/test_grid.py

@@ -1450,3 +1450,64 @@ def test_geo_dataframe(structured_grid, vertex_grid, unstructured_grid):
                     raise AssertionError(
                         f"Cell vertices incorrect for node={node}"
                     )
+
+
+def test_unstructured_iverts_cleanup():
+    grid = GridCases.structured_small()
+
+    # begin building unstructured grid information
+    top = grid.top.ravel()
+    botm = grid.botm[0].ravel()
+    idomain = np.ones(botm.shape, dtype=int)
+
+    # build iac and ja
+    neighbors = grid.neighbors(method="rook", reset=True)
+    iac, ja = [], []
+    for cell, neigh in neighbors.items():
+        iac.append(len(neigh) + 1)
+        ja.extend(
+            [
+                cell,
+            ]
+            + neigh
+        )
+
+    # build iverts and verts without using shared vertices
+    verts, iverts = [], []
+    xverts, yverts = grid.cross_section_vertices
+    ivt = 0
+    for cid, xvs in enumerate(xverts):
+        yvs = yverts[cid]
+        ivts = []
+        for ix, vert in enumerate(xvs[:-1]):
+            ivts.append(ivt)
+            verts.append([ivt, vert, yvs[ix]])
+            ivt += 1
+
+        ivts.append(ivts[0])
+        iverts.append(ivts)
+
+    ugrid = UnstructuredGrid(
+        vertices=verts,
+        iverts=iverts,
+        xcenters=grid.xcellcenters.ravel(),
+        ycenters=grid.ycellcenters.ravel(),
+        iac=iac,
+        ja=ja,
+        top=top,
+        botm=botm,
+        idomain=idomain,
+    )
+
+    if ugrid.nvert != (grid.ncpl * 4):
+        raise AssertionError(
+            "UnstructuredGrid is being built incorrectly for test case"
+        )
+
+    cleaned_vert_num = (grid.nrow + 1) * (grid.ncol + 1)
+    clean_ugrid = ugrid.clean_iverts()
+
+    if clean_ugrid.nvert != cleaned_vert_num:
+        raise AssertionError(
+            "Improper number of vertices for cleaned 'shared' iverts"
+        )

autotest/test_model_splitter.py

@@ -26,7 +26,7 @@ def test_structured_model_splitter(function_tmpdir):
     for row in range(modelgrid.nrow):
         if row != 0 and row % 2 == 0:
             ncol += 1
-        array[row, ncol:] = 2
+        array[row, ncol:] = 100
 
     mfsplit = Mf6Splitter(sim)
     new_sim = mfsplit.split_model(array)
@@ -37,13 +37,13 @@ def test_structured_model_splitter(function_tmpdir):
 
     original_heads = gwf.output.head().get_alldata()[-1]
 
-    ml0 = new_sim.get_model("freyberg_1")
-    ml1 = new_sim.get_model("freyberg_2")
+    ml0 = new_sim.get_model("freyberg_001")
+    ml1 = new_sim.get_model("freyberg_100")
 
     heads0 = ml0.output.head().get_alldata()[-1]
     heads1 = ml1.output.head().get_alldata()[-1]
 
-    new_heads = mfsplit.reconstruct_array({1: heads0, 2: heads1})
+    new_heads = mfsplit.reconstruct_array({1: heads0, 100: heads1})
 
     err_msg = "Heads from original and split models do not match"
     np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg)
@@ -691,3 +691,171 @@ def test_idomain_none(function_tmpdir):
 
     err_msg = "Heads from original and split models do not match"
     np.testing.assert_allclose(new_head, head, atol=1e-07, err_msg=err_msg)
+
+
+@requires_exe("mf6")
+def test_unstructured_complex_disu(function_tmpdir):
+    sim_path = function_tmpdir
+    split_sim_path = sim_path / "model_split"
+
+    # build the simulation structure
+    sim = flopy.mf6.MFSimulation(sim_ws=sim_path)
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    mname = "disu_model"
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=mname)
+
+    # start structured and then create a USG from it
+    nlay = 1
+    nrow = 10
+    ncol = 10
+    delc = np.ones((nrow,))
+    delr = np.ones((ncol,))
+    top = np.ones((nrow, ncol))
+    botm = np.zeros((nlay, nrow, ncol))
+    idomain = np.ones(botm.shape, dtype=int)
+    idomain[0, 1, 4] = 0
+    idomain[0, 8, 5] = 0
+
+    grid = flopy.discretization.StructuredGrid(
+        delc=delc, delr=delr, top=top, botm=botm, idomain=idomain
+    )
+
+    # build the USG connection information
+    neighbors = grid.neighbors(method="rook", reset=True)
+    iac, ja, ihc, cl12, hwva, angldegx = [], [], [], [], [], []
+    for cell, neigh in neighbors.items():
+        iac.append(len(neigh) + 1)
+        ihc.extend(
+            [
+                1,
+            ]
+            * (len(neigh) + 1)
+        )
+        ja.extend(
+            [
+                cell,
+            ]
+            + neigh
+        )
+        cl12.extend(
+            [
+                0,
+            ]
+            + [
+                1,
+            ]
+            * len(neigh)
+        )
+        hwva.extend(
+            [
+                0,
+            ]
+            + [
+                1,
+            ]
+            * len(neigh)
+        )
+        adx = [
+            0,
+        ]
+        for n in neigh:
+            ev = cell - n
+            if ev == -1 * ncol:
+                adx.append(270)
+            elif ev == ncol:
+                adx.append(90)
+            elif ev == -1:
+                adx.append(0)
+            else:
+                adx.append(180)
+        angldegx.extend(adx)
+
+    # build iverts and verts. Do not use shared iverts and mess with verts a
+    #   tiny bit
+    verts, cell2d = [], []
+    xverts, yverts = grid.cross_section_vertices
+    xcenters = grid.xcellcenters.ravel()
+    ycenters = grid.ycellcenters.ravel()
+    ivert = 0
+    for cell_num, xvs in enumerate(xverts):
+        if (cell_num - 3) % 10 == 0:
+            xvs[2] -= 0.001
+            xvs[3] -= 0.001
+        yvs = yverts[cell_num]
+
+        c2drec = [cell_num, xcenters[cell_num], ycenters[cell_num], len(xvs)]
+        for ix, vert in enumerate(xvs[:-1]):
+            c2drec.append(ivert)
+            verts.append([ivert, vert, yvs[ix]])
+            ivert += 1
+
+        c2drec.append(c2drec[4])
+        cell2d.append(c2drec)
+
+    nodes = len(cell2d)
+    nja = len(ja)
+    nvert = len(verts)
+
+    dis = flopy.mf6.ModflowGwfdisu(
+        gwf,
+        nodes=nodes,
+        nja=nja,
+        nvert=nvert,
+        top=np.ravel(grid.top),
+        bot=np.ravel(grid.botm),
+        area=np.ones((nodes,)),
+        idomain=grid.idomain.ravel(),
+        iac=iac,
+        ja=ja,
+        ihc=ihc,
+        cl12=cl12,
+        hwva=hwva,
+        angldegx=angldegx,
+        vertices=verts,
+        cell2d=cell2d,
+    )
+
+    # build npf, ic, CHD, OC package
+    npf = flopy.mf6.ModflowGwfnpf(gwf)
+    ic = flopy.mf6.ModflowGwfic(gwf)
+
+    spd = []
+    for i in range(nrow):
+        spd.append((0 + (i * 10), 0.9))
+        spd.append((9 + (i * 10), 0.5))
+
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=spd)
+
+    spd = {0: [("HEAD", "LAST")]}
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf, head_filerecord=f"{mname}.hds", saverecord=spd
+    )
+
+    sim.write_simulation()
+    sim.run_simulation()
+
+    heads = gwf.output.head().get_alldata()[-1]
+
+    array = np.zeros((nrow, ncol), dtype=int)
+    array[:, 5:] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array)
+
+    new_sim.set_sim_path(split_sim_path)
+    new_sim.write_simulation()
+    new_sim.run_simulation()
+
+    gwf0 = new_sim.get_model(f"{mname}_0")
+    gwf1 = new_sim.get_model(f"{mname}_1")
+
+    heads0 = gwf0.output.head().get_alldata()[-1]
+    heads1 = gwf1.output.head().get_alldata()[-1]
+
+    new_heads = mfsplit.reconstruct_array({0: heads0, 1: heads1})
+
+    diff = np.abs(heads - new_heads)
+    if np.max(diff) > 1e-07:
+        raise AssertionError("Reconstructed head results outside of tolerance")

cliff.toml

@@ -20,13 +20,16 @@ trim = true
 conventional_commits = true
 filter_unconventional = true
 split_commits = false
+commit_preprocessors = [
+    { pattern = "^ *", replace = ""}
+]
 commit_parsers = [
     { message = "^[fF]eat", group = "feat"},
     { message = "^[fF]ix", group = "fix"},
     { message = "^[bB]ug", group = "fix"},
     { message = "^[pP]erf", group = "perf"},
     { message = "^[rR]efactor", group = "refactor"},
-    { message = "^[uU]pdate", group = "refactor"},
+    { message = "^[uU]pdate.*", group = "refactor"},
     { message = "^[dD]oc.*", group = "docs", skip = true},
     { message = "^[bB]inder", group = "docs", skip = true},
     { message = "^[nN]otebook.*", group = "docs", skip = true},

docs/PyPI_release.md

@@ -30,4 +30,4 @@ How to Cite
 
 *Software/Code citation for FloPy:*
 
-[Bakker, Mark, Post, Vincent, Hughes, J. D., Langevin, C. D., White, J. T., Leaf, A. T., Paulinski, S. R., Bellino, J. C., Morway, E. D., Toews, M. W., Larsen, J. D., Fienen, M. N., Starn, J. J., Brakenhoff, D. A., and Bonelli, W. P., 2024, FloPy v3.8.0: U.S. Geological Survey Software Release, 08 August 2024, https://doi.org/10.5066/F7BK19FH](https://doi.org/10.5066/F7BK19FH)
+[Bakker, Mark, Post, Vincent, Hughes, J. D., Langevin, C. D., White, J. T., Leaf, A. T., Paulinski, S. R., Bellino, J. C., Morway, E. D., Toews, M. W., Larsen, J. D., Fienen, M. N., Starn, J. J., Brakenhoff, D. A., and Bonelli, W. P., 2024, FloPy v3.8.1: U.S. Geological Survey Software Release, 05 September 2024, https://doi.org/10.5066/F7BK19FH](https://doi.org/10.5066/F7BK19FH)

flopy/discretization/grid.py

@@ -595,7 +595,7 @@ def xyzvertices(self):
     def cross_section_vertices(self):
         return self.xyzvertices[0], self.xyzvertices[1]
 
-    def geo_dataframe(self, polys):
+    def geo_dataframe(self, features, featuretype="Polygon"):
         """
         Method returns a geopandas GeoDataFrame of the Grid
 
@@ -606,7 +606,7 @@ def geo_dataframe(self, polys):
         from ..utils.geospatial_utils import GeoSpatialCollection
 
         gc = GeoSpatialCollection(
-            polys, shapetype=["Polygon" for _ in range(len(polys))]
+            features, shapetype=[featuretype for _ in range(len(features))]
         )
         gdf = gc.geo_dataframe
         if self.crs is not None:

flopy/discretization/structuredgrid.py

@@ -192,14 +192,18 @@ def __init__(
         if top is not None:
             assert self.__nrow * self.__ncol == len(np.ravel(top))
         if botm is not None:
-            assert self.__nrow * self.__ncol == len(np.ravel(botm[0]))
-            if nlay is not None:
-                self.__nlay = nlay
-            else:
-                if laycbd is not None:
-                    self.__nlay = len(botm) - np.count_nonzero(laycbd)
+            if botm.ndim == 3:
+                assert self.__nrow * self.__ncol == len(np.ravel(botm[0]))
+                if nlay is not None:
+                    self.__nlay = nlay
                 else:
-                    self.__nlay = len(botm)
+                    if laycbd is not None:
+                        self.__nlay = len(botm) - np.count_nonzero(laycbd)
+                    else:
+                        self.__nlay = len(botm)
+            elif botm.ndim == 2:
+                assert botm.shape == (self.__nrow, self.__ncol)
+                self.__nlay = 1
         else:
             self.__nlay = nlay
         if laycbd is not None: