fix(utils/check): file check for mfusg unstructured models (#1145)

* mfusg file checking did not previously work
* add get_neighbors function clause for mfusg unstructured grids (in check class, but not mf6check)
* add modflow/mfdisu/_get_neighboring_nodes function
* fix pakbase/_check_flowp hk and vka for unstructured cases
* modify modflow/mfbas/check function to cater for mfusg unstructured grids, in addition to
  structured modflow grids
* modify t016_test.py to cater for modflow/mfdisu/_get_neighboring_nodes function
* pakbase sidestep np jagged array deprecation warning
* add mfusg load tests to t016

Co-authored-by: Mike Taves <mwtoews@gmail.com>

Resolves #1072

autotest/t016_test.py

@@ -46,11 +46,11 @@ def test_usg_disu_load():
     for (key1, value1), (key2, value2) in zip(
         disu2.__dict__.items(), disu.__dict__.items()
     ):
-        if isinstance(value1, flopy.utils.Util2d) or isinstance(
-            value1, flopy.utils.Util3d
-        ):
+        if isinstance(value1, (flopy.utils.Util2d, flopy.utils.Util3d)):
             assert np.array_equal(value1.array, value2.array)
-        else:
+        elif isinstance(value1, list): #this is for the jagged _get_neighbours list
+            assert np.all([np.all(v1 == v2) for v1,v2 in zip(value1,value2)])
+        elif not isinstance(value1, flopy.utils.reference.TemporalReference):
             assert value1 == value2
 
     return
@@ -136,8 +136,67 @@ def test_usg_model():
             success, buff = mf.run_model()
             assert success
 
+def test_usg_load_01B():
+    print("testing 1-layer unstructured mfusg model loading: 01A_nestedgrid_nognc.nam")
+    pthusgtest = os.path.join(
+        "..", "examples", "data", "mfusg_test", "01A_nestedgrid_nognc"
+    )
+    fname = os.path.join(pthusgtest, "flow.nam")
+    assert os.path.isfile(fname), "nam file not found {}".format(fname)
+
+    # Create the model
+    m = flopy.modflow.Modflow(modelname="usgload_1b", verbose=True)
+
+    # Load the model, with checking
+    m = m.load(fname, check=True)
+
+    # assert disu, lpf, bas packages have been loaded
+    msg = "flopy failed on loading mfusg disu package"
+    assert isinstance(m.disu, flopy.modflow.mfdisu.ModflowDisU), msg
+    msg = "flopy failed on loading mfusg lpf package"
+    assert isinstance(m.lpf, flopy.modflow.mflpf.ModflowLpf), msg
+    msg = "flopy failed on loading mfusg bas package"
+    assert isinstance(m.bas6, flopy.modflow.mfbas.ModflowBas), msg
+    msg = "flopy failed on loading mfusg oc package"
+    assert isinstance(m.oc, flopy.modflow.mfoc.ModflowOc), msg
+    msg = "flopy failed on loading mfusg sms package"
+    assert isinstance(m.sms, flopy.modflow.mfsms.ModflowSms), msg
+
+def test_usg_load_45usg():
+    print("testing 3-layer unstructured mfusg model loading: 45usg.nam")
+    pthusgtest = os.path.join(
+        "..", "examples", "data", "mfusg_test", "45usg"
+    )
+    fname = os.path.join(pthusgtest, "45usg.nam")
+    assert os.path.isfile(fname), "nam file not found {}".format(fname)
+
+    # Create the model
+    m = flopy.modflow.Modflow(modelname="usgload_1b", verbose=True)
+
+    # Load the model, with checking.
+    # RCH not loaded as mfusg rch and evt loading needs work (TO DO)
+    m = m.load(fname, check=True, load_only=["DISU","BAS6","LPF","SMS","OC","DRN","WEL"])
+
+    # assert disu, lpf, bas packages have been loaded
+    msg = "flopy failed on loading mfusg disu package"
+    assert isinstance(m.disu, flopy.modflow.mfdisu.ModflowDisU), msg
+    msg = "flopy failed on loading mfusg lpf package"
+    assert isinstance(m.lpf, flopy.modflow.mflpf.ModflowLpf), msg
+    msg = "flopy failed on loading mfusg bas package"
+    assert isinstance(m.bas6, flopy.modflow.mfbas.ModflowBas), msg
+    msg = "flopy failed on loading mfusg oc package"
+    assert isinstance(m.oc, flopy.modflow.mfoc.ModflowOc), msg
+    msg = "flopy failed on loading mfusg sms package"
+    assert isinstance(m.sms, flopy.modflow.mfsms.ModflowSms), msg
+    msg = "flopy failed on loading mfusg drn package"
+    assert isinstance(m.drn, flopy.modflow.mfdrn.ModflowDrn), msg
+    msg = "flopy failed on loading mfusg wel package"
+    assert isinstance(m.wel, flopy.modflow.mfwel.ModflowWel), msg
+
 
 if __name__ == "__main__":
     test_usg_disu_load()
     test_usg_sms_load()
     test_usg_model()
+    test_usg_load_01B()
+    test_usg_load_45usg()

flopy/modflow/mf.py

@@ -11,6 +11,7 @@
 from ..pakbase import Package
 from ..utils import mfreadnam
 from ..discretization.structuredgrid import StructuredGrid
+from ..discretization.unstructuredgrid import UnstructuredGrid
 from ..discretization.grid import Grid
 from flopy.discretization.modeltime import ModelTime
 from .mfpar import ModflowPar
@@ -264,17 +265,21 @@ def __repr__(self):
 
     @property
     def modeltime(self):
+        if self.get_package("disu") is not None:
+            dis = self.disu
+        else:
+            dis = self.dis
         # build model time
         data_frame = {
-            "perlen": self.dis.perlen.array,
-            "nstp": self.dis.nstp.array,
-            "tsmult": self.dis.tsmult.array,
+            "perlen": dis.perlen.array,
+            "nstp": dis.nstp.array,
+            "tsmult": dis.tsmult.array,
         }
         self._model_time = ModelTime(
             data_frame,
-            self.dis.itmuni_dict[self.dis.itmuni],
-            self.dis.start_datetime,
-            self.dis.steady.array,
+            dis.itmuni_dict[dis.itmuni],
+            dis.start_datetime,
+            dis.steady.array,
         )
         return self._model_time
 
@@ -289,11 +294,18 @@ def modelgrid(self):
             ibound = None
 
         if self.get_package("disu") is not None:
-            self._modelgrid = Grid(
-                grid_type="USG-Unstructured",
-                top=self.disu.top,
-                botm=self.disu.bot,
+            # build unstructured grid
+            self._modelgrid = UnstructuredGrid(
+                grid_type="unstructured",
+                vertices=self._modelgrid.vertices,
+                ivert=self._modelgrid.iverts,
+                xcenters=self._modelgrid.xcenters,
+                ycenters=self._modelgrid.ycenters,
+                ncpl=self.disu.nodelay.array,
+                top=self.disu.top.array,
+                botm=self.disu.bot.array,
                 idomain=ibound,
+                lenuni=self.disu.lenuni,
                 proj4=self._modelgrid.proj4,
                 epsg=self._modelgrid.epsg,
                 xoff=self._modelgrid.xoffset,

flopy/modflow/mfbas.py

@@ -12,7 +12,7 @@
 import sys
 import numpy as np
 from ..pakbase import Package
-from ..utils import Util3d, get_neighbors
+from ..utils import Util3d
 
 
 class ModflowBas(Package):
@@ -215,16 +215,17 @@ def check(self, f=None, verbose=True, level=1, checktype=None):
         """
         chk = self._get_check(f, verbose, level, checktype)
 
-        neighbors = get_neighbors(self.ibound.array)
-        neighbors[
-            np.isnan(neighbors)
-        ] = 0  # set neighbors at edges to 0 (inactive)
-        chk.values(
-            self.ibound.array,
-            (self.ibound.array > 0) & np.all(neighbors < 1, axis=0),
-            "isolated cells in ibound array",
-            "Warning",
-        )
+        neighbors = chk.get_neighbors(self.ibound.array)
+        if isinstance(neighbors, np.ndarray):
+            neighbors[
+                np.isnan(neighbors)
+            ] = 0  # set neighbors at edges to 0 (inactive)
+            chk.values(
+                self.ibound.array,
+                (self.ibound.array > 0) & np.all(neighbors < 1, axis=0),
+                "isolated cells in ibound array",
+                "Warning",
+            )
         chk.values(
             self.ibound.array,
             np.isnan(self.ibound.array),

flopy/modflow/mfdisu.py

@@ -8,6 +8,8 @@
 import numpy as np
 from ..pakbase import Package
 from ..utils import Util2d, Util3d, read1d
+from ..utils.reference import TemporalReference
+from ..discretization.unstructuredgrid import UnstructuredGrid
 
 ITMUNI = {"u": 0, "s": 1, "m": 2, "h": 3, "d": 4, "y": 5}
 LENUNI = {"u": 0, "f": 1, "m": 2, "c": 3}
@@ -453,11 +455,29 @@ def __init__(
             5: "years",
         }
 
+        if start_datetime is None:
+            start_datetime = model._start_datetime
+
+        if model.modelgrid is None:
+            model.modelgrid = UnstructuredGrid(
+                ncpl=self.nodelay.array,
+                top=self.top.array,
+                botm=self.bot.array,
+                lenuni=self.lenuni,
+            )
+
+        self.tr = TemporalReference(
+            itmuni=self.itmuni, start_datetime=start_datetime
+        )
+
         self.start_datetime = start_datetime
 
         # calculate layer thicknesses
         self.__calculate_thickness()
 
+        # get neighboring nodes
+        self._get_neighboring_nodes()
+
         # Add package and return
         self.parent.add_package(self)
         return
@@ -528,7 +548,7 @@ def ncpl(self):
         return self.nodes / self.nlay
 
     @classmethod
-    def load(cls, f, model, ext_unit_dict=None, check=False):
+    def load(cls, f, model, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -924,3 +944,23 @@ def _ftype():
     @staticmethod
     def _defaultunit():
         return 11
+
+    def _get_neighboring_nodes(self):
+        """
+        For each node, get node numbers for all neighbors.
+
+        Returns
+        -------
+        Jagged list of numpy arrays for each node.
+        Each array contains base-1 neighboring node indices.
+        """
+        ja = self.ja.array
+        iac_sum = np.cumsum(self.iac.array)
+        ja_slices = np.asarray(
+            [
+                np.s_[iac_sum[i - 1] + 1 : x] if i > 0 else np.s_[1:x]
+                for i, x in enumerate(iac_sum)
+            ]
+        )  # note: this removes the diagonal - neighbors only
+        self._neighboring_nodes = [ja[sl] for sl in ja_slices]
+        return

flopy/pakbase.py

@@ -112,12 +112,16 @@ def _other_xpf_checks(self, chk, active):
         )
 
         # check vkcb if there are any quasi-3D layers
-        if self.parent.dis.laycbd.sum() > 0:
+        if "DIS" in self.parent.get_package_list():
+            dis = self.parent.dis
+        else:
+            dis = self.parent.disu
+        if dis.laycbd.sum() > 0:
             # pad non-quasi-3D layers in vkcb array with ones so
             # they won't fail checker
             vkcb = self.vkcb.array.copy()
             for l in range(self.vkcb.shape[0]):
-                if self.parent.dis.laycbd[l] == 0:
+                if dis.laycbd[l] == 0:
                     # assign 1 instead of zero as default value that
                     # won't violate checker
                     # (allows for same structure as other checks)
@@ -203,11 +207,21 @@ def _get_kparams(self):
             if kp in self.__dict__:
                 kparams[kp] = name
         if "hk" in self.__dict__:
-            hk = self.hk.array.copy()
+            if self.hk.shape[1] == None:
+                hk = np.asarray(
+                    [a.array.flatten() for a in self.hk], dtype=object
+                )
+            else:
+                hk = self.hk.array.copy()
         else:
             hk = self.k.array.copy()
         if "vka" in self.__dict__ and self.layvka.sum() > 0:
-            vka = self.vka.array
+            if self.vka.shape[1] == None:
+                vka = np.asarray(
+                    [a.array.flatten() for a in self.vka], dtype=object
+                )
+            else:
+                vka = self.vka.array
             vka_param = kparams.pop("vka")
         elif "k33" in self.__dict__:
             vka = self.k33.array

flopy/utils/__init__.py

@@ -48,7 +48,7 @@
     SwrListBudget,
     Mf6ListBudget,
 )
-from .check import check, get_neighbors
+from .check import check
 from .utils_def import FlopyBinaryData, totim_to_datetime
 from .flopy_io import read_fixed_var, write_fixed_var
 from .zonbud import (