From 345e0fbfdc007b688982c09bfa82841d56ac5cf7 Mon Sep 17 00:00:00 2001 From: Chris Nicol <37314969+cnicol-gwlogic@users.noreply.github.com> Date: Sat, 10 Jul 2021 06:59:03 +1000 Subject: [PATCH] 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 Resolves #1072 --- autotest/t016_test.py | 67 ++++++++++++++++++++++++-- flopy/modflow/mf.py | 32 +++++++++---- flopy/modflow/mfbas.py | 23 ++++----- flopy/modflow/mfdisu.py | 42 ++++++++++++++++- flopy/pakbase.py | 22 +++++++-- flopy/utils/__init__.py | 2 +- flopy/utils/check.py | 102 ++++++++++++++++++++++++++-------------- 7 files changed, 225 insertions(+), 65 deletions(-) diff --git a/autotest/t016_test.py b/autotest/t016_test.py index a30c4f64e..c5a3627c0 100644 --- a/autotest/t016_test.py +++ b/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() diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py index e16b942a8..19d92beef 100644 --- a/flopy/modflow/mf.py +++ b/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, diff --git a/flopy/modflow/mfbas.py b/flopy/modflow/mfbas.py index 35c30cfae..a5fb35ae7 100644 --- a/flopy/modflow/mfbas.py +++ b/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), diff --git a/flopy/modflow/mfdisu.py b/flopy/modflow/mfdisu.py index 4682c486a..4aaf63a90 100644 --- a/flopy/modflow/mfdisu.py +++ b/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 diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 69c453956..0c7f219a8 100644 --- a/flopy/pakbase.py +++ b/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 diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 26f815ea5..91dcc1134 100644 --- a/flopy/utils/__init__.py +++ b/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 ( diff --git a/flopy/utils/check.py b/flopy/utils/check.py index e594c2380..683f7e434 100644 --- a/flopy/utils/check.py +++ b/flopy/utils/check.py @@ -2,6 +2,7 @@ import numpy as np from numpy.lib import recfunctions from ..utils.recarray_utils import recarray +from ..utils.util_array import Util3d class check: @@ -666,6 +667,73 @@ def _get_dtype(self): ] ) + def get_neighbors(self, a): + """ + For a structured grid, this returns the 6 neighboring values for each + value in a. For an unstructured grid, this returns the n_max neighboring + values, where n_max is the maximum number of nodal connections for any + node within the model; nodes with less than n_max connections are assigned + np.nan for indices above the number of connections for that node. + + Parameters + ---------- + a : 3-D Model array in layer, row, column order array, even for an + unstructured grid; for instance, a Util3d array + (e.g. flopy.modflow.ModflowBas.ibound). + + Returns + ------- + neighbors : 4-D array + Array of neighbors, where axis 0 contains the n neighboring + values for each value in a, and subsequent axes are in layer, row, + column order. "n" is 6 for a structured grid, and "n" is n_max + for an unstructured grid, as described above. Nan is returned for + values at edges. + """ + if self.structured: + nk, ni, nj = a.shape + tmp = np.empty((nk + 2, ni + 2, nj + 2), dtype=float) + tmp[:, :, :] = np.nan + tmp[1:-1, 1:-1, 1:-1] = a[:, :, :] + neighbors = np.vstack( + [ + tmp[0:-2, 1:-1, 1:-1].ravel(), # k-1 + tmp[2:, 1:-1, 1:-1].ravel(), # k+1 + tmp[1:-1, 0:-2, 1:-1].ravel(), # i-1 + tmp[1:-1, 2:, 1:-1].ravel(), # i+1 + tmp[1:-1, 1:-1, :-2].ravel(), # j-1 + tmp[1:-1, 1:-1, 2:].ravel(), + ] + ) # j+1 + return neighbors.reshape(6, nk, ni, nj) + else: + if "DISU" in self.model.get_package_list(): + disu = self.model.disu + neighbors = disu._neighboring_nodes + if isinstance(a, Util3d): + a = a.array + pad_value = int(-1e9) + n_max = ( + np.max(disu.iac.array) - 1 + ) # -1 for self, removed below + arr_neighbors = [ + np.pad( + a[n - 1], + (0, n_max - n.size), + "constant", + constant_values=(0, pad_value), + ) + for n in neighbors + ] + arr_neighbors = np.where( + arr_neighbors == -1e9, np.nan, arr_neighbors + ) + neighbors = arr_neighbors.T + else: + # if no disu, we can't define neighbours for this ugrid + neighbors = None + return neighbors + def _fmt_string_list(array, float_format="{}"): fmt_string = [] @@ -738,40 +806,6 @@ def fields_view(arr, fields): return np.ndarray(arr.shape, dtype2, arr, 0, arr.strides) -def get_neighbors(a): - """ - Returns the 6 neighboring values for each value in a. - - Parameters - ---------- - a : 3-D array - Model array in layer, row, column order. - - Returns - ------- - neighbors : 4-D array - Array of neighbors, where axis 0 contains the 6 neighboring - values for each value in a, and subsequent axes are in layer, row, - column order. - Nan is returned for values at edges. - """ - nk, ni, nj = a.shape - tmp = np.empty((nk + 2, ni + 2, nj + 2), dtype=float) - tmp[:, :, :] = np.nan - tmp[1:-1, 1:-1, 1:-1] = a[:, :, :] - neighbors = np.vstack( - [ - tmp[0:-2, 1:-1, 1:-1].ravel(), # k-1 - tmp[2:, 1:-1, 1:-1].ravel(), # k+1 - tmp[1:-1, 0:-2, 1:-1].ravel(), # i-1 - tmp[1:-1, 2:, 1:-1].ravel(), # i+1 - tmp[1:-1, 1:-1, :-2].ravel(), # j-1 - tmp[1:-1, 1:-1, 2:].ravel(), - ] - ) # j+1 - return neighbors.reshape(6, nk, ni, nj) - - class mf6check(check): def __init__( self,