feat(lgr-disv): add to_disv_gridprops() method to lgr object (modflow…

…py#2271)

* feat(lgr-disv): add new to_disv_gridprops to lgr object

* ruff ruff

* modflow6 lgr tutorial was mistakenly marked as modflow-lgr tutorial

* add test

* add layer information to gridprops and update test

* ruff

* minor tweaks to require idomain the same for all layers

* add deprecation warnings for gridlist_to_disv_gridprops

* ruff

* ruff

.docs/tutorials.rst

@@ -35,6 +35,7 @@ MODFLOW 6
    Notebooks/mf6_output_tutorial01
    Notebooks/mf6_sfr_tutorial01
    Notebooks/mf6_tutorial01
+   Notebooks/mf6_lgr_tutorial01
 
 
 MODFLOW-2005
@@ -50,15 +51,6 @@ MODFLOW-2005
    Notebooks/mf_tutorial02
 
 
-MODFLOW-LGR
------------
-
-.. toctree::
-   :maxdepth: 2
-
-   Notebooks/lgr_tutorial01
-
-
 MODFLOW-NWT
 -----------
 

autotest/test_export.py

@@ -1718,50 +1718,51 @@ def test_vtk_export_disv1_model(function_tmpdir):
         idomain=np.ones((nlay, nrow, ncol)),
     )
 
-    from flopy.utils.cvfdutil import gridlist_to_disv_gridprops
-
-    gridprops = gridlist_to_disv_gridprops([mg])
-    gridprops["top"] = 0
-    gridprops["botm"] = np.zeros((nlay, nrow * ncol), dtype=float) - 1
-    gridprops["nlay"] = nlay
-
-    disv = ModflowGwfdisv(gwf, **gridprops)
-    ic = ModflowGwfic(gwf, strt=10)
-    npf = ModflowGwfnpf(gwf)
-
-    # Export model without specifying packages_names parameter
-    # create the vtk output
-    gwf = sim.get_model()
-    vtkobj = Vtk(gwf, binary=False)
-    vtkobj.add_model(gwf)
-    f = function_tmpdir / "gwf.vtk"
-    vtkobj.write(f)
-
-    # load the output using the vtk standard library
-    gridreader = vtkUnstructuredGridReader()
-    gridreader.SetFileName(str(f))
-    gridreader.Update()
-    grid = gridreader.GetOutput()
-
-    # get the points
-    vtk_points = grid.GetPoints()
-    vtk_points = vtk_points.GetData()
-    vtk_points = vtk_to_numpy(vtk_points)
-
-    # get cell locations (ia format of point to cell relationship)
-    cell_locations = vtk_to_numpy(grid.GetCellLocationsArray())
-    cell_locations_answer = np.array([0, 8, 16, 24, 32, 40, 48, 56, 64])
-    print(f"Found cell locations {cell_locations} in vtk file.")
-    print(f"Expecting cell locations {cell_locations_answer}")
-    errmsg = "vtk cell locations do not match expected result."
-    assert np.allclose(cell_locations, cell_locations_answer), errmsg
-
-    cell_types = vtk_to_numpy(grid.GetCellTypesArray())
-    cell_types_answer = np.array(9 * [42])
-    print(f"Found cell types {cell_types} in vtk file.")
-    print(f"Expecting cell types {cell_types_answer}")
-    errmsg = "vtk cell types do not match expected result."
-    assert np.allclose(cell_types, cell_types_answer), errmsg
+    with pytest.deprecated_call():
+        from flopy.utils.cvfdutil import gridlist_to_disv_gridprops
+
+        gridprops = gridlist_to_disv_gridprops([mg])
+        gridprops["top"] = 0
+        gridprops["botm"] = np.zeros((nlay, nrow * ncol), dtype=float) - 1
+        gridprops["nlay"] = nlay
+
+        disv = ModflowGwfdisv(gwf, **gridprops)
+        ic = ModflowGwfic(gwf, strt=10)
+        npf = ModflowGwfnpf(gwf)
+
+        # Export model without specifying packages_names parameter
+        # create the vtk output
+        gwf = sim.get_model()
+        vtkobj = Vtk(gwf, binary=False)
+        vtkobj.add_model(gwf)
+        f = function_tmpdir / "gwf.vtk"
+        vtkobj.write(f)
+
+        # load the output using the vtk standard library
+        gridreader = vtkUnstructuredGridReader()
+        gridreader.SetFileName(str(f))
+        gridreader.Update()
+        grid = gridreader.GetOutput()
+
+        # get the points
+        vtk_points = grid.GetPoints()
+        vtk_points = vtk_points.GetData()
+        vtk_points = vtk_to_numpy(vtk_points)
+
+        # get cell locations (ia format of point to cell relationship)
+        cell_locations = vtk_to_numpy(grid.GetCellLocationsArray())
+        cell_locations_answer = np.array([0, 8, 16, 24, 32, 40, 48, 56, 64])
+        print(f"Found cell locations {cell_locations} in vtk file.")
+        print(f"Expecting cell locations {cell_locations_answer}")
+        errmsg = "vtk cell locations do not match expected result."
+        assert np.allclose(cell_locations, cell_locations_answer), errmsg
+
+        cell_types = vtk_to_numpy(grid.GetCellTypesArray())
+        cell_types_answer = np.array(9 * [42])
+        print(f"Found cell types {cell_types} in vtk file.")
+        print(f"Expecting cell types {cell_types_answer}")
+        errmsg = "vtk cell types do not match expected result."
+        assert np.allclose(cell_types, cell_types_answer), errmsg
 
 
 @pytest.mark.mf6

autotest/test_grid.py

@@ -943,25 +943,27 @@ def test_tocvfd3():
         yoff=200,
         idomain=idomain,
     )
-    gridprops = gridlist_to_disv_gridprops([sg1, sg2])
-    assert "ncpl" in gridprops
-    assert "nvert" in gridprops
-    assert "vertices" in gridprops
-    assert "cell2d" in gridprops
-
-    ncpl = gridprops["ncpl"]
-    nvert = gridprops["nvert"]
-    vertices = gridprops["vertices"]
-    cell2d = gridprops["cell2d"]
-    assert ncpl == 121
-    assert nvert == 148
-    assert len(vertices) == nvert
-    assert len(cell2d) == 121
-
-    # spot check information for cell 28 (zero based)
-    answer = [28, 250.0, 150.0, 7, 38, 142, 143, 45, 46, 44, 38]
-    for i, j in zip(cell2d[28], answer):
-        assert i == j, f"{i} not equal {j}"
+
+    with pytest.deprecated_call():
+        gridprops = gridlist_to_disv_gridprops([sg1, sg2])
+        assert "ncpl" in gridprops
+        assert "nvert" in gridprops
+        assert "vertices" in gridprops
+        assert "cell2d" in gridprops
+
+        ncpl = gridprops["ncpl"]
+        nvert = gridprops["nvert"]
+        vertices = gridprops["vertices"]
+        cell2d = gridprops["cell2d"]
+        assert ncpl == 121
+        assert nvert == 148
+        assert len(vertices) == nvert
+        assert len(cell2d) == 121
+
+        # spot check information for cell 28 (zero based)
+        answer = [28, 250.0, 150.0, 7, 38, 142, 143, 45, 46, 44, 38]
+        for i, j in zip(cell2d[28], answer):
+            assert i == j, f"{i} not equal {j}"
 
 
 @requires_pkg("shapely")

autotest/test_lgrutil.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from flopy.utils.lgrutil import Lgr
+from flopy.utils.lgrutil import Lgr, LgrToDisv
 
 
 def test_lgrutil():
@@ -155,3 +155,75 @@ def test_lgrutil2():
     ]
     assert np.allclose(lgr.delr, answer), f"{lgr.delr} /= {answer}"
     assert np.allclose(lgr.delc, answer), f"{lgr.delc} /= {answer}"
+
+
+def test_lgrutil3():
+    # Define parent grid information
+    xoffp = 0.0
+    yoffp = 0.0
+    nlayp = 3
+    nrowp = 3
+    ncolp = 3
+
+    dx = 100.0
+    dy = 100.0
+    dz = 10.0
+    delrp = dx * np.ones(ncolp)
+    delcp = dy * np.ones(nrowp)
+    topp = dz * np.ones((nrowp, ncolp), dtype=float)
+    botmp = np.empty((nlayp, nrowp, ncolp), dtype=float)
+    for k in range(nlayp):
+        botmp[k] = -(k + 1) * dz
+    idomainp = np.ones((nlayp, nrowp, ncolp), dtype=int)
+    idomainp[:, nrowp // 2, ncolp // 2] = 0
+    ncpp = 3
+    ncppl = nlayp * [1]
+    lgr = Lgr(
+        nlayp,
+        nrowp,
+        ncolp,
+        delrp,
+        delcp,
+        topp,
+        botmp,
+        idomainp,
+        ncpp=ncpp,
+        ncppl=ncppl,
+        xllp=xoffp,
+        yllp=yoffp,
+    )
+
+    # check to make sure gridprops is accessible from lgr
+    gridprops = lgr.to_disv_gridprops()
+    assert "ncpl" in gridprops
+    assert "nvert" in gridprops
+    assert "vertices" in gridprops
+    assert "nlay" in gridprops
+    assert "top" in gridprops
+    assert "botm" in gridprops
+    assert gridprops["ncpl"] == 17
+    assert gridprops["nvert"] == 32
+    assert gridprops["nlay"] == 3
+
+    # test the lgr to disv class
+    lgrtodisv = LgrToDisv(lgr)
+
+    # test guts of LgrToDisv to make sure hanging vertices added correctly
+    assert lgrtodisv.right_face_hanging[(1, 0)] == [0, 4, 8, 12]
+    assert lgrtodisv.left_face_hanging[(1, 2)] == [3, 7, 11, 15]
+    assert lgrtodisv.back_face_hanging[(2, 1)] == [12, 13, 14, 15]
+    assert lgrtodisv.front_face_hanging[(0, 1)] == [0, 1, 2, 3]
+
+    assert lgrtodisv.iverts[1] == [1, 2, 6, 18, 17, 5]
+    assert lgrtodisv.iverts[3] == [4, 5, 20, 24, 9, 8]
+    assert lgrtodisv.iverts[4] == [6, 7, 11, 10, 27, 23]
+    assert lgrtodisv.iverts[6] == [9, 29, 30, 10, 14, 13]
+
+    assert np.allclose(gridprops["top"], dz * np.ones((17,)))
+
+    assert gridprops["botm"].shape == (3, 17)
+    b = np.empty((3, 17))
+    b[0] = -dz
+    b[1] = -2 * dz
+    b[2] = -3 * dz
+    assert np.allclose(gridprops["botm"], b)

flopy/utils/cvfdutil.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import pandas as pd
 
@@ -390,6 +392,10 @@ def gridlist_to_disv_gridprops(gridlist):
     be numbered according to consecutive numbering of active cells in the
     grid list.
 
+    This function is deprecated in 3.8 and will be removed in 3.9.  Use the
+    functionality in flopy.utils.cvfdutil.Lgr() to create a DISV mesh for a
+    nested grid.
+
     Parameters
     ----------
     gridlist : list
@@ -403,6 +409,13 @@ def gridlist_to_disv_gridprops(gridlist):
         modflow6 disv package.
 
     """
+    warnings.warn(
+        "the gridlist_to_disv_gridprops function is deprecated and will be "
+        "removed in version 3.9. Use flopy.utils.cvfdutil.Lgr() instead, which "
+        "allows a nested grid to be created and exported to a DISV mesh.",
+        PendingDeprecationWarning,
+    )
+
     verts, iverts = gridlist_to_verts(gridlist)
     gridprops = get_disv_gridprops(verts, iverts)
     return gridprops