update(optimize_splitting_mask): add active_only flag

autotest/test_model_splitter.py

@@ -221,7 +221,10 @@ def test_metis_splitting_with_lak_sfr(function_tmpdir):
 @requires_exe("mf6")
 @requires_pkg("pymetis")
 @requires_pkg("h5py")
+# @requires_pkg("sklearn")
 def test_save_load_node_mapping_structured(function_tmpdir):
+    import pymetis
+
     sim_path = get_example_data_path() / "mf6-freyberg"
     new_sim_path = function_tmpdir / "mf6-freyberg/split_model"
     hdf_file = new_sim_path / "node_map.hdf5"
@@ -235,7 +238,9 @@ def test_save_load_node_mapping_structured(function_tmpdir):
     original_heads = sim.get_model().output.head().get_alldata()[-1]
 
     mfsplit = Mf6Splitter(sim)
-    array = mfsplit.optimize_splitting_mask(nparts=nparts)
+    array = mfsplit.optimize_splitting_mask(
+        nparts=nparts, active_only=True, options=pymetis.Options(seed=42, contig=1)
+    )
     new_sim = mfsplit.split_model(array)
     new_sim.set_sim_path(new_sim_path)
     new_sim.write_simulation()

etc/environment.yml

@@ -64,3 +64,4 @@ dependencies:
   - vtk
   - xmipy
   - h5py
+  - scikit-learn

flopy/discretization/grid.py

@@ -71,7 +71,7 @@ class Grid:
         The value can be anything accepted by
         :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
         such as an authority string (eg "EPSG:26916") or a WKT string.
-    prjfile : str or PathLike, optional if `crs` is specified
+    prjfile : str or pathlike, optional if `crs` is specified
         ESRI-style projection file with well-known text defining the CRS
         for the model grid (must be projected; geographic CRS are not supported).
     xoff : float
@@ -648,7 +648,7 @@ def _set_neighbors(self, reset=False, method="rook"):
         ----------
         reset : bool
             flag to recalculate neighbors
-        method: str
+        method : str
             "rook" for shared edges and "queen" for shared vertex
 
         Returns
@@ -722,6 +722,7 @@ def neighbors(self, node=None, **kwargs):
         """
         method = kwargs.pop("method", None)
         reset = kwargs.pop("reset", False)
+
         if method is None:
             self._set_neighbors(reset=reset)
         else:
@@ -744,6 +745,32 @@ def neighbors(self, node=None, **kwargs):
 
         return self._neighbors
 
+    def get_shared_edge(self, node0, node1, iverts=True):
+        """
+        Method to get the shared iverts or vertices between two cells. The shared
+        edge defines the shared cell face.
+
+        Parameters
+        ----------
+        node0 : int
+        node1 : int
+        iverts : bool
+            boolean flag to return shared iverts (True) or shared vertices (False)
+
+        Returns
+        -------
+            tuple : iverts or vertices that define the shared face
+        """
+        iv0 = set(self.iverts[node0])
+        iv1 = set(self.iverts[node1])
+        shared = tuple(iv0 & iv1)
+
+        if iverts:
+            return tuple(shared)
+        else:
+            verts = [self.verts[shared[0]], self.verts[shared[1]]]
+            return tuple(verts)
+
     def remove_confining_beds(self, array):
         """
         Method to remove confining bed layers from an array

flopy/discretization/structuredgrid.py

@@ -590,15 +590,15 @@ def is_regular_z(self):
                 0, 0, 0
             ]
             failed = np.abs(rel_diff_thick0) > rel_tol
-            is_regular_z = np.count_nonzero(failed) == 0
+            _is_regular_z = np.count_nonzero(failed) == 0
             for k in range(1, self.nlay):
                 rel_diff_zk = (self.delz[k, :, :] - self.delz[0, :, :]) / self.delz[
                     0, :, :
                 ]
                 failed = np.abs(rel_diff_zk) > rel_tol
-                is_regular_z = is_regular_z and np.count_nonzero(failed) == 0
+                _is_regular_z = _is_regular_z and np.count_nonzero(failed) == 0
 
-            self._cache_dict[cache_index] = CachedData(is_regular_z)
+            self._cache_dict[cache_index] = CachedData(_is_regular_z)
         if self._copy_cache:
             return self._cache_dict[cache_index].data
         else:
@@ -801,7 +801,7 @@ def convert_grid(self, factor):
             raise AssertionError("Grid is not complete and cannot be converted")
 
     ###############
-    # Methods
+    ### Methods ###
     ###############
     def neighbors(self, *args, **kwargs):
         """
@@ -835,6 +835,7 @@ def neighbors(self, *args, **kwargs):
         """
         nn = None
         as_nodes = kwargs.pop("as_nodes", False)
+        memhog = kwargs.pop("fast", False)
 
         if kwargs:
             if "node" in kwargs:
@@ -864,12 +865,105 @@ def neighbors(self, *args, **kwargs):
         else:
             as_nodes = True
 
-        neighbors = super().neighbors(nn, **kwargs)
-        if not as_nodes:
-            neighbors = self.get_lrc(neighbors)
+        if memhog:
+            neighbors = self._fast_neighbors(**kwargs)
+        else:
+            neighbors = super().neighbors(nn, **kwargs)
+            if not as_nodes:
+                neighbors = self.get_lrc(neighbors)
 
         return neighbors
 
+    def _fast_neighbors(self, **kwargs):
+        """
+        Memory intensive and not elegent in any sense, but very fast method to get
+        neighbors for Structured Grids
+
+        Returns
+        -------
+        dict
+        """
+        from collections import defaultdict
+
+        if self._neighbors is None or kwargs.pop("reset", False):
+            nrow = self.nrow
+            ncol = self.ncol
+
+            ncpl = nrow * ncol
+
+            arr = np.arange(ncpl, dtype=np.int32).reshape((nrow, ncol))
+
+            # general case
+            arr2 = arr[1:-1, 1:-1].ravel()
+            neighs2 = np.empty((4, (nrow - 2) * (ncol - 2)), dtype=np.int32)
+            neighs2[0] = arr2 - ncol  # u
+            neighs2[1] = arr2 + 1  # r
+            neighs2[2] = arr2 + ncol  # d
+            neighs2[3] = arr2 - 1  # l
+            neighs2 = neighs2.T.tolist()
+            neighbors = {v: neighs2[ix] for ix, v in enumerate(arr2)}
+
+            # ecase 1
+            arr2 = arr[0, 1:-1].ravel()
+            neighs2 = np.empty((3, ncol - 2), dtype=np.int32)
+            # no up
+            neighs2[0] = arr2 + 1  # r
+            neighs2[1] = arr2 + ncol  # d
+            neighs2[2] = arr2 - 1  # l
+            neighs2 = neighs2.T.tolist()
+
+            for ix, v in enumerate(arr2):
+                neighbors[v] = neighs2[ix]
+
+            # ecase 2
+            arr2 = arr[-1, 1:-1].ravel()
+            neighs2 = np.empty((3, ncol - 2), dtype=np.int32)
+            neighs2[0] = arr2 - ncol  # u
+            neighs2[1] = arr2 + 1  # r
+            # no down
+            neighs2[2] = arr2 - 1  # l
+            neighs2 = neighs2.T.tolist()
+
+            for ix, v in enumerate(arr2):
+                neighbors[v] = neighs2[ix]
+
+            # ecase 3
+            arr2 = arr[1:-1, 0].ravel()
+            neighs2 = np.empty((3, nrow - 2), dtype=np.int32)
+            neighs2[0] = arr2 - ncol  # u
+            neighs2[1] = arr2 + 1  # r
+            neighs2[2] = arr2 + ncol  # d
+            # no left
+            neighs2 = neighs2.T.tolist()
+
+            for ix, v in enumerate(arr2):
+                neighbors[v] = neighs2[ix]
+
+            # ecase 4
+            arr2 = arr[1:-1, -1].ravel()
+            neighs2 = np.empty((3, nrow - 2), dtype=np.int32)
+            neighs2[0] = arr2 - ncol  # u
+            # no right
+            neighs2[1] = arr2 + ncol  # d
+            neighs2[2] = arr2 - 1  # l
+
+            neighs2 = neighs2.T.tolist()
+
+            for ix, v in enumerate(arr2):
+                neighbors[v] = neighs2[ix]
+
+            del arr
+
+            # corners
+            neighbors[0] = [1, ncol]  # r, d
+            neighbors[ncol - 1] = [(ncol - 1) + ncol, ncol - 2]  # d, l
+            neighbors[ncpl - ncol] = [(ncpl - ncol) - ncol, (ncpl - ncol) + 1]  # u, r
+            neighbors[ncpl - 1] = [(ncpl - 1) - ncol, ncpl - 2]  # u, l
+
+            self._neighbors = neighbors
+
+        return self._neighbors
+
     def intersect(self, x, y, z=None, local=False, forgive=False):
         """
         Get the row and column of a point with coordinates x and y
@@ -1682,11 +1776,15 @@ def _set_structured_iverts(self):
         pair for each vertex
 
         """
-        iverts = []
-        for i in range(self.nrow):
-            for j in range(self.ncol):
-                iverts.append(self._build_structured_iverts(i, j))
-        self._iverts = iverts
+        rowarr = np.repeat(np.arange(self.nrow, dtype=int), self.ncol)
+        colarr = np.tile(np.arange(self.ncol, dtype=int), self.nrow)
+
+        iverts = np.empty((4, self.ncpl), dtype=int)
+        iverts[0] = rowarr * (self.ncol + 1) + colarr
+        iverts[1] = rowarr * (self.ncol + 1) + colarr + 1
+        iverts[2] = (rowarr + 1) * (self.ncol + 1) + colarr + 1
+        iverts[3] = (rowarr + 1) * (self.ncol + 1) + colarr
+        self._iverts = iverts.T.tolist()
         return
 
     def _build_structured_iverts(self, i, j):