From 7378eaeb4e27560e44c5d0bda73ab20b20af5bb4 Mon Sep 17 00:00:00 2001 From: Chris Nicol Date: Tue, 6 Jul 2021 08:16:30 +1000 Subject: [PATCH 1/5] fix(utils/check) file check for mfusg unstructured models 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 Resolves #1072 --- autotest/t016_test.py | 8 ++-- flopy/modflow/mf.py | 32 +++++++++----- flopy/modflow/mfbas.py | 35 +++++++-------- flopy/modflow/mfdisu.py | 41 +++++++++++++++++- flopy/pakbase.py | 20 ++++++--- flopy/utils/__init__.py | 2 +- flopy/utils/check.py | 95 ++++++++++++++++++++++++++--------------- 7 files changed, 160 insertions(+), 73 deletions(-) diff --git a/autotest/t016_test.py b/autotest/t016_test.py index a30c4f64ef..bf9876e4ad 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 diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py index e16b942a8c..19d92beef6 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 35c30cfae8..a96fd1ba58 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,22 +215,23 @@ 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", - ) - chk.values( - self.ibound.array, - np.isnan(self.ibound.array), - error_name="Not a number", - error_type="Error", - ) + 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), + error_name="Not a number", + error_type="Error", + ) chk.summarize() return chk diff --git a/flopy/modflow/mfdisu.py b/flopy/modflow/mfdisu.py index 4682c486a4..b2037d6da7 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,28 @@ 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 +547,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 +943,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 69c453956d..d0f4ecb16d 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,17 @@ 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]) + 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]) + else: + vka = self.vka.array vka_param = kparams.pop("vka") elif "k33" in self.__dict__: vka = self.k33.array @@ -215,7 +225,7 @@ def _get_kparams(self): else: vka = None if vka is not None: - vka = vka.copy() + vka = vka.copy() return kparams, hk, vka, vka_param def _check_flowp(self, f=None, verbose=True, level=1, checktype=None): diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 26f815ea55..91dcc1134b 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 e594c23807..f36521efa6 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,65 @@ 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 = [] @@ -737,41 +797,6 @@ def fields_view(arr, fields): dtype2 = np.dtype({name: arr.dtype.fields[name] for name in 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, From 49023e88dad5a72387c980f616eba734132e9eca Mon Sep 17 00:00:00 2001 From: Chris Nicol Date: Tue, 6 Jul 2021 08:18:49 +1000 Subject: [PATCH 2/5] run black --- flopy/modflow/mfdisu.py | 3 ++- flopy/pakbase.py | 2 +- flopy/utils/check.py | 21 +++++++++++++++------ 3 files changed, 18 insertions(+), 8 deletions(-) diff --git a/flopy/modflow/mfdisu.py b/flopy/modflow/mfdisu.py index b2037d6da7..4aaf63a90d 100644 --- a/flopy/modflow/mfdisu.py +++ b/flopy/modflow/mfdisu.py @@ -463,7 +463,8 @@ def __init__( ncpl=self.nodelay.array, top=self.top.array, botm=self.bot.array, - lenuni=self.lenuni) + lenuni=self.lenuni, + ) self.tr = TemporalReference( itmuni=self.itmuni, start_datetime=start_datetime diff --git a/flopy/pakbase.py b/flopy/pakbase.py index d0f4ecb16d..de1b71bffc 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -225,7 +225,7 @@ def _get_kparams(self): else: vka = None if vka is not None: - vka = vka.copy() + vka = vka.copy() return kparams, hk, vka, vka_param def _check_flowp(self, f=None, verbose=True, level=1, checktype=None): diff --git a/flopy/utils/check.py b/flopy/utils/check.py index f36521efa6..683f7e434e 100644 --- a/flopy/utils/check.py +++ b/flopy/utils/check.py @@ -680,7 +680,7 @@ def get_neighbors(self, a): 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 @@ -713,20 +713,28 @@ def get_neighbors(self, a): if isinstance(a, Util3d): a = a.array pad_value = int(-1e9) - n_max = np.max(disu.iac.array) - 1 # -1 for self, removed below + 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] + 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) + 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 = [] for field in array.dtype.descr: @@ -797,6 +805,7 @@ def fields_view(arr, fields): dtype2 = np.dtype({name: arr.dtype.fields[name] for name in fields}) return np.ndarray(arr.shape, dtype2, arr, 0, arr.strides) + class mf6check(check): def __init__( self, From e6687168f10e8232c8df5d69d50d7dba74567ef2 Mon Sep 17 00:00:00 2001 From: Chris Nicol <37314969+cnicol-gwlogic@users.noreply.github.com> Date: Tue, 6 Jul 2021 09:14:58 +1000 Subject: [PATCH 3/5] Update autotest/t016_test.py Co-authored-by: Mike Taves --- autotest/t016_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/autotest/t016_test.py b/autotest/t016_test.py index bf9876e4ad..eec1a2d145 100644 --- a/autotest/t016_test.py +++ b/autotest/t016_test.py @@ -46,7 +46,7 @@ def test_usg_disu_load(): for (key1, value1), (key2, value2) in zip( disu2.__dict__.items(), disu.__dict__.items() ): - if isinstance(value1, (flopy.utils.Util2d, flopy.utils.Util3d): + if isinstance(value1, (flopy.utils.Util2d, flopy.utils.Util3d)): assert np.array_equal(value1.array, value2.array) 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)]) From ce85a010d57b4f6884ac214f8955ee6df1bb417c Mon Sep 17 00:00:00 2001 From: Chris Nicol Date: Tue, 6 Jul 2021 11:04:26 +1000 Subject: [PATCH 4/5] pakbase sidestep np jagged array deprecation warning --- flopy/pakbase.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/flopy/pakbase.py b/flopy/pakbase.py index de1b71bffc..0c7f219a80 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -208,14 +208,18 @@ def _get_kparams(self): kparams[kp] = name if "hk" in self.__dict__: if self.hk.shape[1] == None: - hk = np.asarray([a.array.flatten() for a in self.hk]) + 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: if self.vka.shape[1] == None: - vka = np.asarray([a.array.flatten() for a in self.vka]) + vka = np.asarray( + [a.array.flatten() for a in self.vka], dtype=object + ) else: vka = self.vka.array vka_param = kparams.pop("vka") From 7a6dc0b55a2bf224519084ff677f59e8da33d3f3 Mon Sep 17 00:00:00 2001 From: Chris Nicol Date: Wed, 7 Jul 2021 09:17:08 +1000 Subject: [PATCH 5/5] * add mfusg load tests to t016 --- autotest/t016_test.py | 59 ++++++++++++++++++++++++++++++++++++++++++ flopy/modflow/mfbas.py | 12 ++++----- 2 files changed, 65 insertions(+), 6 deletions(-) diff --git a/autotest/t016_test.py b/autotest/t016_test.py index eec1a2d145..c5a3627c0f 100644 --- a/autotest/t016_test.py +++ b/autotest/t016_test.py @@ -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/mfbas.py b/flopy/modflow/mfbas.py index a96fd1ba58..a5fb35ae7d 100644 --- a/flopy/modflow/mfbas.py +++ b/flopy/modflow/mfbas.py @@ -226,12 +226,12 @@ def check(self, f=None, verbose=True, level=1, checktype=None): "isolated cells in ibound array", "Warning", ) - chk.values( - self.ibound.array, - np.isnan(self.ibound.array), - error_name="Not a number", - error_type="Error", - ) + chk.values( + self.ibound.array, + np.isnan(self.ibound.array), + error_name="Not a number", + error_type="Error", + ) chk.summarize() return chk