feat(dis2d): introduce limited support for a 2D structured grid (for …

…overland flow) (#2131)

* This first PR doesn't really do anything, except define model shapes and cellids for dis2d
* No testing yet
* A modelgrid for this discretization type has not been addressed yet
* Unclear what type of plotting and export to support

flopy/mf6/coordinates/modeldimensions.py

@@ -405,6 +405,10 @@ def _create_model_grid(self, grid_type):
             self._model_grid = ModelGrid(
                 self.model_name, self.simulation_data, DiscretizationType.DISL
             )
+        elif grid_type == DiscretizationType.DIS2D:
+            self._model_grid = ModelGrid(
+                self.model_name, self.simulation_data, DiscretizationType.DIS2D
+            )
         else:
             self._model_grid = ModelGrid(
                 self.model_name,

flopy/mf6/coordinates/modelgrid.py

@@ -457,7 +457,13 @@ def get_grid_type(simulation_data, model_name):
             is not None
         ):
             return DiscretizationType.DISL
-
+        elif (
+            package_recarray.search_data(
+                f"dis2d{structure.get_version_string()}", 0
+            )
+            is not None
+        ):
+            return DiscretizationType.DIS2D
         return DiscretizationType.UNDEFINED
 
     def get_idomain(self):
@@ -477,6 +483,10 @@ def get_idomain(self):
             return self._simulation_data.mfdata[
                 (self._model_name, "disu", "griddata", "idomain")
             ].get_data()
+        elif self._grid_type == DiscretizationType.DIS2D:
+            return self._simulation_data.mfdata[
+                (self._model_name, "dis2d", "griddata", "idomain")
+            ].get_data()
         except_str = (
             "ERROR: Grid type {} for model {} not " "recognized.".format(
                 self._grid_type, self._model_name
@@ -529,6 +539,8 @@ def get_horizontal_cross_section_dim_arrays(self):
     def get_model_dim(self):
         if self.grid_type() == DiscretizationType.DIS:
             return [self.num_layers(), self.num_rows(), self.num_columns()]
+        elif self.grid_type() == DiscretizationType.DIS2D:
+            return [self.num_rows(), self.num_columns()]
         elif self.grid_type() == DiscretizationType.DISV:
             return [self.num_layers(), self.num_cells_per_layer()]
         elif (
@@ -544,6 +556,11 @@ def get_model_dim_arrays(self):
                 np.arange(1, self.num_rows() + 1, 1, np.int32),
                 np.arange(1, self.num_columns() + 1, 1, np.int32),
             ]
+        elif self.grid_type() == DiscretizationType.DIS2D:
+            return [
+                np.arange(1, self.num_rows() + 1, 1, np.int32),
+                np.arange(1, self.num_columns() + 1, 1, np.int32),
+            ]
         elif self.grid_type() == DiscretizationType.DISV:
             return [
                 np.arange(1, self.num_layers() + 1, 1, np.int32),
@@ -584,6 +601,8 @@ def get_horizontal_cross_section_dim_names(self):
     def get_model_dim_names(self):
         if self.grid_type() == DiscretizationType.DIS:
             return ["layer", "row", "column"]
+        elif self.grid_type() == DiscretizationType.DIS2D:
+            return ["row", "column"]
         elif self.grid_type() == DiscretizationType.DISV:
             return ["layer", "layer_cell_num"]
         elif (
@@ -596,38 +615,48 @@ def get_num_spatial_coordinates(self):
         grid_type = self.grid_type()
         if grid_type == DiscretizationType.DIS:
             return 3
+        elif grid_type == DiscretizationType.DIS2D:
+            return 2
         elif grid_type == DiscretizationType.DISV:
             return 2
-        elif (
-            grid_type == DiscretizationType.DISU
-            or grid_type == DiscretizationType.DISL
-        ):
+        elif grid_type == DiscretizationType.DISU:
+            return 1
+        elif grid_type == DiscretizationType.DISL:
             return 1
+        return 0
 
     def num_rows(self):
-        if self.grid_type() != DiscretizationType.DIS:
+        if self.grid_type() not in [
+            DiscretizationType.DIS,
+            DiscretizationType.DIS2D,
+        ]:
             except_str = (
                 'ERROR: Model "{}" does not have rows.  Can not '
                 "return number of rows.".format(self._model_name)
             )
             print(except_str)
             raise MFGridException(except_str)
 
+        distype = self.grid_type().name.lower()  # dis or dis2d
         return self._simulation_data.mfdata[
-            (self._model_name, "dis", "dimensions", "nrow")
+            (self._model_name, distype, "dimensions", "nrow")
         ].get_data()
 
     def num_columns(self):
-        if self.grid_type() != DiscretizationType.DIS:
+        if self.grid_type() not in [
+            DiscretizationType.DIS,
+            DiscretizationType.DIS2D,
+        ]:
             except_str = (
                 'ERROR: Model "{}" does not have columns.  Can not '
                 "return number of columns.".format(self._model_name)
             )
             print(except_str)
             raise MFGridException(except_str)
 
+        distype = self.grid_type().name.lower()  # dis or dis2d
         return self._simulation_data.mfdata[
-            (self._model_name, "dis", "dimensions", "ncol")
+            (self._model_name, distype, "dimensions", "ncol")
         ].get_data()
 
     def num_connections(self):
@@ -668,12 +697,15 @@ def num_layers(self):
         elif (
             self.grid_type() == DiscretizationType.DISU
             or self.grid_type() == DiscretizationType.DISL
+            or self.grid_type() == DiscretizationType.DIS2D
         ):
             return None
 
     def num_cells(self):
         if self.grid_type() == DiscretizationType.DIS:
             return self.num_rows() * self.num_columns() * self.num_layers()
+        elif self.grid_type() == DiscretizationType.DIS2D:
+            return self.num_rows() * self.num_columns()
         elif self.grid_type() == DiscretizationType.DISV:
             return self.num_layers() * self.num_cells_per_layer()
         elif self.grid_type() == DiscretizationType.DISU:
@@ -693,6 +725,11 @@ def get_all_model_cells(self):
                     for column in range(0, self.num_columns()):
                         model_cells.append((layer + 1, row + 1, column + 1))
             return model_cells
+        elif self.grid_type() == DiscretizationType.DIS2D:
+            for row in range(0, self.num_rows()):
+                for column in range(0, self.num_columns()):
+                    model_cells.append((layer + 1, row + 1, column + 1))
+            return model_cells
         elif self.grid_type() == DiscretizationType.DISV:
             for layer in range(0, self.num_layers()):
                 for layer_cellid in range(0, self.num_rows()):

flopy/mf6/utils/mfenums.py

@@ -11,3 +11,4 @@ class DiscretizationType(Enum):
     DISV = 2
     DISU = 3
     DISL = 4
+    DIS2D = 5