From 088f147f1d761875f3a1ffabdda82312c002305a Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Thu, 16 Jan 2020 07:26:35 -0800 Subject: [PATCH] feat(mf6 checker): input data check for mf6 (#779) * feat(input data check for mf6): refactored mf2005 input data check to support both mf2005 and mf6 * feat(input data check for mf6): refactor of mf2005 code to work with mf6 * feat(input data check for mf6): autotest updates * feat(input data check for mf6) * fix(codacy fixes) * fix(check too complex) * fix(codacy): simplified check code --- autotest/t024_test.py | 1 + autotest/t504_test.py | 27 +- autotest/t505_test.py | 46 ++- flopy/discretization/structuredgrid.py | 10 + flopy/discretization/unstructuredgrid.py | 7 +- flopy/export/netcdf.py | 3 +- flopy/mbase.py | 157 +++++--- flopy/mf6/data/mfdataarray.py | 25 +- flopy/mf6/data/mfdatalist.py | 15 + flopy/mf6/data/mfdatascalar.py | 15 + flopy/mf6/data/mfdatastorage.py | 51 ++- flopy/mf6/data/mfdatautil.py | 2 +- flopy/mf6/data/mffileaccess.py | 4 +- flopy/mf6/mfmodel.py | 40 +- flopy/mf6/mfpackage.py | 39 +- flopy/mf6/modflow/mfsimulation.py | 71 +++- flopy/modflow/mf.py | 6 +- flopy/modflow/mfbas.py | 6 +- flopy/modflow/mfdis.py | 6 +- flopy/modflow/mfmnw2.py | 14 +- flopy/modflow/mfmnwi.py | 6 +- flopy/modflow/mfoc.py | 5 +- flopy/modflow/mfrch.py | 7 +- flopy/modflow/mfriv.py | 9 +- flopy/modflow/mfsfr2.py | 2 +- flopy/modpath/mpsim.py | 6 +- flopy/mt3d/mt.py | 8 +- flopy/pakbase.py | 492 +++++++++++++---------- flopy/seawat/swt.py | 6 +- flopy/utils/check.py | 254 +++++++++--- 30 files changed, 920 insertions(+), 420 deletions(-) diff --git a/autotest/t024_test.py b/autotest/t024_test.py index 1e355592b..2aa488d19 100644 --- a/autotest/t024_test.py +++ b/autotest/t024_test.py @@ -40,6 +40,7 @@ def test_bcs_check(): ibound[1, 1, 1] = 1 # fully isolated cell ibound[0:2, 4, 4] = 1 # cell connected vertically to one other cell bas = flopy.modflow.ModflowBas(mf, ibound=ibound) + mf._mg_resync = True chk = bas.check() assert chk.summary_array['desc'][0] == 'isolated cells in ibound array' assert chk.summary_array.i[0] == 1 and chk.summary_array.i[0] == 1 and \ diff --git a/autotest/t504_test.py b/autotest/t504_test.py index 177966f61..4561d1c35 100644 --- a/autotest/t504_test.py +++ b/autotest/t504_test.py @@ -49,7 +49,7 @@ def test001a_tharmonic(): # load simulation sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, - verbosity_level=0) + verbosity_level=0, verify_data=True) sim.simulation_data.mfpath.set_sim_path(run_folder) # write simulation to new location @@ -150,7 +150,8 @@ def test003_gwfs_disv(): array_util = PyListUtil() # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data=True) # make temp folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -235,7 +236,7 @@ def test005_advgw_tidal(): # load simulation sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, - verbosity_level=2) + verbosity_level=2, verify_data=True) # test obs/ts package interface model = sim.get_model(model_name) @@ -321,7 +322,8 @@ def test006_gwf3(): array_util = PyListUtil() # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data=True) model = sim.get_model() disu = model.get_package('disu') @@ -460,7 +462,8 @@ def test045_lake1ss_table(): expected_head_file_b = os.path.join(expected_output_folder, 'lakeex1b_adj.hds') # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data=True) # make temp folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -532,7 +535,7 @@ def test006_2models_mvr(): expected_head_file_bb = os.path.join(expected_output_folder, 'model2_adj.hds') # load simulation - sim = MFSimulation.load(sim_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(sim_name, 'mf6', exe_name, pth, verify_data=True) # make temp folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -685,7 +688,8 @@ def test001e_uzf_3lay(): expected_head_file_b = os.path.join(expected_output_folder, 'test001e_UZF_3lay_adj.hds') # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data=True) # make temp folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -766,7 +770,8 @@ def test045_lake2tr(): 'lakeex2a_adj.hds') # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data=True) # write simulation to new location sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -824,7 +829,8 @@ def test036_twrihfb(): expected_head_file_b = os.path.join(expected_output_folder, 'twrihfb2015_output_adj.hds') # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data = True) # make temp folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -890,7 +896,8 @@ def test027_timeseriestest(): expected_head_file_b = os.path.join(expected_output_folder, 'timeseriestest_adj.hds') # load simulation - sim = MFSimulation.load(model_name, 'mf6', exe_name, pth) + sim = MFSimulation.load(model_name, 'mf6', exe_name, pth, + verify_data=True) # make temp folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) diff --git a/autotest/t505_test.py b/autotest/t505_test.py index 104f4fcfd..c1dfe4745 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -74,8 +74,8 @@ def np001(): exe_name=exe_name, sim_ws=run_folder, continue_=True, memory_print_option='summary') name = test_sim.name_file - assert name.continue_.get_data() == True - assert name.nocheck.get_data() == None + assert name.continue_.get_data() + assert name.nocheck.get_data() is None assert name.memory_print_option.get_data() == 'summary' kwargs = {} @@ -289,8 +289,25 @@ def np001(): stress_period_data=well_spd) except MFDataException: error_occurred = True - assert error_occurred + + # test error checking + drn_package = ModflowGwfdrn(model, print_input=True, print_flows=True, + save_flows=True, maxbound=1, + timeseries=[(0.0, 60.0), (100000.0, 60.0)], + stress_period_data=[((100, 0, 0), np.nan, + 'drn_1'), ((0, 0, 0), + 10.0)]) + npf_package = ModflowGwfnpf(model, save_flows=True, + alternative_cell_averaging='logarithmic', + icelltype=1, k=100001.0, k33=1e-12) + chk = sim.check() + summary = '.'.join(chk[0].summary_array.desc) + assert 'drn_1 package: invalid BC index' in summary + assert 'npf package: vertical hydraulic conductivity values below ' \ + 'checker threshold of 1e-11' in summary + assert 'npf package: horizontal hydraulic conductivity values above ' \ + 'checker threshold of 100000.0' in summary return @@ -391,8 +408,8 @@ def np002(): sim.run_simulation() sim2 = MFSimulation.load(sim_ws=run_folder) - model = sim2.get_model(model_name) - npf_package = model.get_package('npf') + model_ = sim2.get_model(model_name) + npf_package = model_.get_package('npf') k = npf_package.k.array # get expected results @@ -426,6 +443,20 @@ def np002(): # clean up sim.delete_output_files() + # test error checking + sto_package = ModflowGwfsto(model, save_flows=True, iconvert=1, + ss=0.00000001, sy=0.6) + chd_package = ModflowGwfchd(model, print_input=True, print_flows=True, + maxbound=1, stress_period_data=[((0, 0, 0), + np.nan)]) + chk = sim.check() + summary = '.'.join(chk[0].summary_array.desc) + assert 'sto package: specific storage values below ' \ + 'checker threshold of 1e-06' in summary + assert 'sto package: specific yield values above ' \ + 'checker threshold of 0.5' in summary + assert 'Not a number' in summary + return @@ -954,6 +985,11 @@ def test005_advgw_tidal(): # clean up sim.delete_output_files() + # check packages + chk = sim.check() + summary = '.'.join(chk[0].summary_array.desc) + assert summary == '' + return diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 524c98e05..beaa663d7 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -55,9 +55,15 @@ def __init__(self, delc=None, delr=None, top=None, botm=None, idomain=None, assert self.__nrow * self.__ncol == len(np.ravel(top)) if botm is not None: assert self.__nrow * self.__ncol == len(np.ravel(botm[0])) + if nlay is not None and nlay < len(botm): + self.__nlay_nocbd = nlay + else: + self.__nlay_nocbd = len(botm) + self.__nlay = len(botm) else: self.__nlay = nlay + self.__nlay_nocbd = nlay #################### # Properties @@ -75,6 +81,10 @@ def is_complete(self): return True return False + @property + def nlay_nocbd(self): + return self.__nlay_nocbd + @property def nlay(self): return self.__nlay diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index 0885d2172..3e71061e3 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -96,11 +96,16 @@ def nnodes(self): @property def ncpl(self): if self._ncpl is None: - return len(self._iverts) + if self._iverts is None: + return None + else: + return len(self._iverts) return self._ncpl @property def shape(self): + if self.ncpl is None: + return self.nnodes if isinstance(self.ncpl, (list, np.ndarray)): return self.nlay, self.ncpl[0] else: diff --git a/flopy/export/netcdf.py b/flopy/export/netcdf.py index d0df25d58..9fdf22a23 100644 --- a/flopy/export/netcdf.py +++ b/flopy/export/netcdf.py @@ -197,7 +197,8 @@ def __init__(self, output_filename, model, time_values=None, self.z_positive = z_positive if self.grid_units is None: self.grid_units = "unspecified" - assert self.grid_units in ["feet", "meters", "unspecified"], \ + assert self.grid_units in ["feet", "meters", "unspecified", + "undefined"], \ "unsupported length units: " + self.grid_units self.time_units = self.model_time.time_units diff --git a/flopy/mbase.py b/flopy/mbase.py index ea154fef7..7271f7855 100644 --- a/flopy/mbase.py +++ b/flopy/mbase.py @@ -158,6 +158,105 @@ def verbose(self): 'must define verbose in child ' 'class to use this base class') + @abc.abstractmethod + def check(self, f=None, verbose=True, level=1): + raise NotImplementedError( + 'must define check in child ' + 'class to use this base class') + + def get_package_list(self, ftype=None): + """ + Get a list of all the package names. + + Parameters + ---------- + ftype : str + Type of package, 'RIV', 'LPF', etc. + + Returns + ------- + val : list of strings + Can be used to see what packages are in the model, and can then + be used with get_package to pull out individual packages. + + """ + val = [] + for pp in (self.packagelist): + if ftype is None: + val.append(pp.name[0].upper()) + elif pp.package_type.lower() == ftype: + val.append(pp.name[0].upper()) + return val + + def _check(self, chk, level=1): + """ + Check model data for common errors. + + Parameters + ---------- + f : str or file handle + String defining file name or file handle for summary file + of check method output. If a string is passed a file handle + is created. If f is None, check method does not write + results to a summary file. (default is None) + verbose : bool + Boolean flag used to determine if check method results are + written to the screen + level : int + Check method analysis level. If level=0, summary checks are + performed. If level=1, full checks are performed. + summarize : bool + Boolean flag used to determine if summary of results is written + to the screen + + Returns + ------- + None + + Examples + -------- + + >>> import flopy + >>> m = flopy.modflow.Modflow.load('model.nam') + >>> m.check() + """ + + # check instance for model-level check + results = {} + + for p in self.packagelist: + if chk.package_check_levels.get(p.name[0].lower(), 0) <= level: + results[p.name[0]] = p.check(f=None, verbose=False, + level=level - 1, + checktype=chk.__class__) + + # model level checks + # solver check + if self.version in chk.solver_packages.keys(): + solvers = set(chk.solver_packages[self.version]).intersection( + set(self.get_package_list())) + if not solvers: + chk._add_to_summary('Error', desc='\r No solver package', + package='model') + elif len(list(solvers)) > 1: + for s in solvers: + chk._add_to_summary('Error', + desc='\r Multiple solver packages', + package=s) + else: + chk.passed.append('Compatible solver package') + + # add package check results to model level check summary + for r in results.values(): + if r is not None and r.summary_array is not None: # currently SFR doesn't have one + chk.summary_array = np.append(chk.summary_array, + r.summary_array).view( + np.recarray) + chk.passed += ['{} package: {}'.format(r.package.name[0], psd) + for psd in r.passed] + chk.summarize() + return chk + class BaseModel(ModelInterface): """ @@ -1003,30 +1102,6 @@ def get_package(self, name): return pp return None - def get_package_list(self, ftype=None): - """ - Get a list of all the package names. - - Parameters - ---------- - ftype : str - Type of package, 'RIV', 'LPF', etc. - - Returns - ------- - val : list of strings - Can be used to see what packages are in the model, and can then - be used with get_package to pull out individual packages. - - """ - val = [] - for pp in (self.packagelist): - if ftype is None: - val.append(pp.name[0].upper()) - elif pp.package_type.lower() == ftype: - val.append(pp.name[0].upper()) - return val - def set_version(self, version): self.version = version.lower() @@ -1367,29 +1442,6 @@ def check(self, f=None, verbose=True, level=1): # check instance for model-level check chk = utils.check(self, f=f, verbose=verbose, level=level) - results = {} - - for p in self.packagelist: - if chk.package_check_levels.get(p.name[0].lower(), 0) <= level: - results[p.name[0]] = p.check(f=None, verbose=False, - level=level - 1) - - # model level checks - # solver check - if self.version in chk.solver_packages.keys(): - solvers = set(chk.solver_packages[self.version]).intersection( - set(self.get_package_list())) - if not solvers: - chk._add_to_summary('Error', desc='\r No solver package', - package='model') - elif len(list(solvers)) > 1: - for s in solvers: - chk._add_to_summary('Error', - desc='\r Multiple solver packages', - package=s) - else: - chk.passed.append('Compatible solver package') - # check for unit number conflicts package_units = {} duplicate_units = {} @@ -1408,16 +1460,7 @@ def check(self, f=None, verbose=True, level=1): else: chk.passed.append('Unit number conflicts') - # add package check results to model level check summary - for k, r in results.items(): - if r is not None and r.summary_array is not None: # currently SFR doesn't have one - chk.summary_array = np.append(chk.summary_array, - r.summary_array).view( - np.recarray) - chk.passed += ['{} package: {}'.format(r.package.name[0], psd) - for psd in r.passed] - chk.summarize() - return chk + return self._check(chk, level) def plot(self, SelPackList=None, **kwargs): """ diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index 4074fd9f2..b507ec399 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -32,6 +32,17 @@ class MFArray(MFMultiDimVar): dimensions : MFDataDimensions dimension information related to the model, package, and array + Attributes + ---------- + data_type : DataType + type of data stored in the scalar + plotable : bool + if the scalar is plotable + dtype : numpy.dtype + the scalar's numpy data type + data : variable + calls get_data with default parameters + Methods ------- new_simulation : (sim_data : MFSimulationData) @@ -369,8 +380,6 @@ def make_layered(self): def store_as_external_file(self, external_file_path, multiplier=None, layer=None, binary=False): - if multiplier is None: - multiplier = [1.0] if isinstance(layer, int): layer = (layer,) storage = self._get_storage_obj() @@ -378,6 +387,8 @@ def store_as_external_file(self, external_file_path, multiplier=None, self._set_storage_obj(self._new_storage(False, True)) storage = self._get_storage_obj() ds_index = self._resolve_layer_index(layer) + if multiplier is None: + multiplier = [storage.get_default_mult()] try: # move data to file @@ -428,6 +439,10 @@ def has_data(self, layer=None): value_, traceback_, None, self._simulation_data.debug, ex) + @property + def data(self): + return self.get_data() + def get_data(self, layer=None, apply_mult=False, **kwargs): if self._get_storage_obj() is None: self._data_storage = self._new_storage(False) @@ -454,11 +469,11 @@ def get_data(self, layer=None, apply_mult=False, **kwargs): return None def set_data(self, data, multiplier=None, layer=None): - if multiplier is None: - multiplier = [1.0] self._resync() if self._get_storage_obj() is None: self._data_storage = self._new_storage(False) + if multiplier is None: + multiplier = [self._get_storage_obj().get_default_mult()] if isinstance(layer, int): layer = (layer,) if isinstance(data, str): @@ -992,8 +1007,6 @@ def get_data(self, layer=None, apply_mult=True, **kwargs): return None def set_data(self, data, multiplier=None, layer=None, key=None): - if multiplier is None: - multiplier = [1.0] if isinstance(data, dict) or isinstance(data, OrderedDict): # each item in the dictionary is a list for one stress period # the dictionary key is the stress period the list is for diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index 4871a02e0..4798a19d7 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -34,6 +34,17 @@ class MFList(mfdata.MFMultiDimVar, DataListInterface): dimensions : MFDataDimensions dimension information related to the model, package, and array + Attributes + ---------- + data_type : DataType + type of data stored in the scalar + plotable : bool + if the scalar is plotable + dtype : numpy.dtype + the scalar's numpy data type + data : variable + calls get_data with default parameters + Methods ------- new_simulation : (sim_data : MFSimulationData) @@ -929,6 +940,10 @@ def add_transient_key(self, transient_key): self._data_storage[transient_key] = \ super(MFTransientList, self)._new_storage(stress_period) + @property + def data(self): + return self.get_data() + def get_data(self, key=None, apply_mult=False, **kwargs): if self._data_storage is not None and len(self._data_storage) > 0: if key is None: diff --git a/flopy/mf6/data/mfdatascalar.py b/flopy/mf6/data/mfdatascalar.py index e493f17fe..8eca36c80 100644 --- a/flopy/mf6/data/mfdatascalar.py +++ b/flopy/mf6/data/mfdatascalar.py @@ -30,6 +30,17 @@ class MFScalar(mfdata.MFData): dimensions : MFDataDimensions dimension information related to the model, package, and array + Attributes + ---------- + data_type : DataType + type of data stored in the scalar + plotable : bool + if the scalar is plotable + dtype : numpy.dtype + the scalar's numpy data type + data : variable + calls get_data with default parameters + Methods ------- has_data : () : bool @@ -107,6 +118,10 @@ def has_data(self): value_, traceback_, None, self._simulation_data.debug, ex) + @property + def data(self): + return self.get_data() + def get_data(self, apply_mult=False, **kwargs): try: return self._get_storage_obj().get_data(apply_mult=apply_mult) diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py index cd1eed72a..ff0e7ea4e 100644 --- a/flopy/mf6/data/mfdatastorage.py +++ b/flopy/mf6/data/mfdatastorage.py @@ -86,14 +86,19 @@ class LayerStorage(object): """ def __init__(self, data_storage, lay_indexes, - data_storage_type=DataStorageType.internal_array): + data_storage_type=DataStorageType.internal_array, + data_type=None): self._data_storage_parent = data_storage self._lay_indexes = lay_indexes self.internal_data = None self.data_const_value = None self.data_storage_type = data_storage_type + self.data_type = data_type self.fname = None - self.factor = 1.0 + if self.data_type == DatumType.integer: + self.factor = 1 + else: + self.factor = 1.0 self.iprn = None self.binary = False @@ -274,6 +279,11 @@ def __init__(self, sim_data, model_or_sim, data_dimensions, get_file_entry, self._data_storage_type = data_storage_type self._stress_period = stress_period self._data_path = data_path + if not data_structure_type == DataStructureType.recarray: + self._data_type = self.data_dimensions.structure.\ + get_datum_type(return_enum_type=True) + else: + self._data_type = None self.layer_storage = MultiList(shape=layer_shape, callback=self._create_layer) #self.layer_storage = [LayerStorage(self, x, data_storage_type) @@ -287,10 +297,7 @@ def __init__(self, sim_data, model_or_sim, data_dimensions, get_file_entry, if data_structure_type == DataStructureType.recarray: self.build_type_list(resolve_data_shape=False) - self._data_type = None - else: - self._data_type = self.data_dimensions.structure.\ - get_datum_type(return_enum_type=True) + self.layered = layered # initialize comments @@ -304,13 +311,15 @@ def __str__(self): return self.get_data_str(False) def _create_layer(self, indexes): - return LayerStorage(self, indexes, self._data_storage_type) + return LayerStorage(self, indexes, self._data_storage_type, + self._data_type) def flatten(self): self.layered = False storage_type = self.layer_storage.first_item().data_storage_type self.layer_storage = MultiList(mdlist=[LayerStorage(self, 0, - storage_type)]) + storage_type, + self._data_type)]) def make_layered(self): if not self.layered: @@ -881,7 +890,7 @@ def store_internal(self, data, layer=None, const=False, multiplier=None, key=None, autofill=False, print_format=None): if multiplier is None: - multiplier = [1.0] + multiplier = [self.get_default_mult()] if self.data_structure_type == DataStructureType.recarray: if self.layer_storage.first_item().data_storage_type == \ DataStorageType.internal_constant: @@ -1029,7 +1038,7 @@ def store_external(self, file_path, layer=None, multiplier=None, print_format=None, data=None, do_not_verify=False, binary=False): if multiplier is None: - multiplier = [1.0] + multiplier = [self.get_default_mult()] layer_new, multiplier = self._store_prep(layer, multiplier) if data is not None: @@ -1304,10 +1313,13 @@ def _validate_cellid(self, arr_line, data_index): cellid_size = model_grid.get_num_spatial_coordinates() if cellid_size + data_index > len(arr_line): return False - for index in range(data_index, cellid_size + data_index): + for index, \ + dim_size in zip(range(data_index, cellid_size + data_index), + model_grid.get_model_dim()): if not DatumUtil.is_int(arr_line[index]): return False - if int(arr_line[index]) <= 0: + val = int(arr_line[index]) + if val <= 0 or val > dim_size: return False return True @@ -1323,10 +1335,7 @@ def add_data_line_comment(self, comment, line_num): line_num) def process_internal_line(self, arr_line): - if self._data_type == DatumType.integer: - multiplier = 1 - else: - multiplier = 1.0 + multiplier = self.get_default_mult() print_format = None if isinstance(arr_line, list): if len(arr_line) < 2: @@ -1952,6 +1961,12 @@ def build_type_list(self, data_set=None, data=None, (data_item.name, data_type)) return self._recarray_type_list + def get_default_mult(self): + if self._data_type == DatumType.integer: + return 1 + else: + return 1.0 + @staticmethod def _calc_data_size(data, count_to=None, current_length=None): if current_length is None: @@ -2020,12 +2035,12 @@ def _store_prep(self, layer, multiplier): mult_ml = MultiList(multiplier) if not mult_ml.in_shape(layer): if multiplier[0] is None: - multiplier = 1.0 + multiplier = self.get_default_mult() else: multiplier = multiplier[0] else: if mult_ml.first_item() is None: - multiplier = 1.0 + multiplier = self.get_default_mult() else: multiplier = mult_ml.first_item() diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 70336aee6..a6eaa4b48 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -601,7 +601,7 @@ def empty(self, model, maxbound=None, aux_vars=None, boundnames=False, for aux_var in aux_vars: type_list.append((aux_var, object)) if boundnames: - type_list.append(('boundnames', object)) + type_list.append(('boundname', object)) if timeseries: # fix type list to make all types objects diff --git a/flopy/mf6/data/mffileaccess.py b/flopy/mf6/data/mffileaccess.py index 5025f6c67..c697f4bc9 100644 --- a/flopy/mf6/data/mffileaccess.py +++ b/flopy/mf6/data/mffileaccess.py @@ -514,7 +514,7 @@ def _load_layer(self, layer, layer_size, storage, arr_line, file_handle, try: storage.store_internal([convert_data( arr_line[1], self._data_dimensions, self.structure.type, - di_struct)], layer, const=True, multiplier=[1.0]) + di_struct)], layer, const=True) except Exception as ex: type_, value_, traceback_ = sys.exc_info() raise MFDataException(self.structure.get_model(), @@ -867,7 +867,7 @@ def read_list_data_from_file(self, file_handle, storage, current_key, convert_data(arr_line[1], self._data_dimensions, struct.data_item_structures[1].type, struct.data_item_structures[0]), - 0, const=True, multiplier=[1.0]) + 0, const=True) else: data_rec = storage._build_recarray(arr_line[1], None, True) diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index f46e8c1be..55965aab3 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -18,6 +18,7 @@ from ..mbase import ModelInterface from .utils.mfenums import DiscretizationType from .data import mfstructure +from ..utils.check import mf6check class MFModel(PackageContainer, ModelInterface): @@ -365,7 +366,8 @@ def modelgrid(self): botm=dis.bot.array, idomain=idomain, lenuni=dis.length_units.array, proj4=self._modelgrid.proj4, epsg=self._modelgrid.epsg, xoff=self._modelgrid.xoffset, - yoff=self._modelgrid.yoffset, angrot=self._modelgrid.angrot) + yoff=self._modelgrid.yoffset, angrot=self._modelgrid.angrot, + nodes=dis.nodes.get_data()) else: return self._modelgrid @@ -464,6 +466,40 @@ def verbose(self): def verbose(self, verbose): self._verbose = verbose + def check(self, f=None, verbose=True, level=1): + """ + Check model data for common errors. + + Parameters + ---------- + f : str or file handle + String defining file name or file handle for summary file + of check method output. If a string is passed a file handle + is created. If f is None, check method does not write + results to a summary file. (default is None) + verbose : bool + Boolean flag used to determine if check method results are + written to the screen + level : int + Check method analysis level. If level=0, summary checks are + performed. If level=1, full checks are performed. + + Returns + ------- + None + + Examples + -------- + + >>> import flopy + >>> m = flopy.modflow.Modflow.load('model.nam') + >>> m.check() + """ + # check instance for model-level check + chk = mf6check(self, f=f, verbose=verbose, level=level) + + return self._check(chk, level) + @classmethod def load_base(cls, simulation, structure, modelname='NewModel', model_nam_file='modflowtest.nam', mtype='gwf', version='mf6', @@ -800,7 +836,7 @@ def remove_package(self, package_name): if not isinstance(packages, list): packages = [packages] for package in packages: - if package._model_or_sim.name != self.name: + if package.model_or_sim.name != self.name: except_text = 'Package can not be removed from model {} ' \ 'since it is ' \ 'not part of ' diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index c357e8286..ded786646 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -16,6 +16,7 @@ from .coordinates import modeldimensions from ..pakbase import PackageInterface from .data.mfdatautil import MFComment +from ..utils.check import mf6check class MFBlockHeader(object): @@ -1121,7 +1122,7 @@ class MFPackage(PackageContainer, PackageInterface): def __init__(self, model_or_sim, package_type, filename=None, pname=None, loading_package=False, parent_file=None): - self._model_or_sim = model_or_sim + self.model_or_sim = model_or_sim self._data_list = [] self._package_type = package_type if model_or_sim.type == 'Model' and package_type.lower() != 'nam': @@ -1166,7 +1167,7 @@ def __init__(self, model_or_sim, package_type, filename=None, pname=None, if filename is None: self._filename = MFFileMgmt.string_to_file_path('{}.{}'.format( - self._model_or_sim.name, package_type)) + self.model_or_sim.name, package_type)) else: if not isinstance(filename, str): message = 'Invalid fname parameter. Expecting type str. ' \ @@ -1264,7 +1265,7 @@ def parent(self, parent): @property def plotable(self): - if self._model_or_sim.type == "Simulation": + if self.model_or_sim.type == "Simulation": return False else: return True @@ -1274,15 +1275,29 @@ def data_list(self): # return [data_object, data_object, ...] return self._data_list + def check(self, f=None, verbose=True, level=1, checktype=None): + if checktype is None: + checktype = mf6check + return super(MFPackage, self).check(f, verbose, level, checktype) + + def _get_nan_exclusion_list(self): + excl_list = [] + if hasattr(self, 'stress_period_data'): + spd_struct = self.stress_period_data.structure + for item_struct in spd_struct.data_item_structures: + if item_struct.optional or item_struct.keystring_dict: + excl_list.append(item_struct.name) + return excl_list + def _get_data_str(self, formal, show_data=True): data_str = 'package_name = {}\nfilename = {}\npackage_type = {}' \ '\nmodel_or_simulation_package = {}' \ '\n{}_name = {}' \ '\n'.format(self._get_pname(), self._filename, self.package_type, - self._model_or_sim.type.lower(), - self._model_or_sim.type.lower(), - self._model_or_sim.name) + self.model_or_sim.type.lower(), + self.model_or_sim.type.lower(), + self.model_or_sim.name) if self.parent_file is not None and formal: data_str = '{}parent_file = ' \ '{}\n\n'.format(data_str, self.parent_file._get_pname()) @@ -1408,16 +1423,16 @@ def _update_size_defs(self): dataset.structure.name)) def remove(self): - self._model_or_sim.remove_package(self) + self.model_or_sim.remove_package(self) def build_child_packages_container(self, pkg_type, filerecord): # get package class package_obj = self.package_factory(pkg_type, - self._model_or_sim.model_type) + self.model_or_sim.model_type) # create child package object child_pkgs_name = 'utl{}packages'.format(pkg_type) child_pkgs_obj = self.package_factory(child_pkgs_name, '') - child_pkgs = child_pkgs_obj(self._model_or_sim, self, pkg_type, + child_pkgs = child_pkgs_obj(self.model_or_sim, self, pkg_type, filerecord, None, package_obj) setattr(self, pkg_type, child_pkgs) self._child_package_groups[pkg_type] = child_pkgs @@ -1431,8 +1446,8 @@ def build_child_package(self, pkg_type, data, parameter_name, filerecord): child_path = package_group._next_default_file_path() # create new empty child package package_obj = self.package_factory(pkg_type, - self._model_or_sim.model_type) - package = package_obj(self._model_or_sim, filename=child_path, + self.model_or_sim.model_type) + package = package_obj(self.model_or_sim, filename=child_path, parent_file=self) assert hasattr(package, parameter_name) @@ -1472,7 +1487,7 @@ def build_mfdata(self, var_name, data=None): self.blocks[block.name] = MFBlock(self._simulation_data, self.dimensions, block, self.path + (key,), - self._model_or_sim, self) + self.model_or_sim, self) dataset_struct = block.data_structures[var_name] var_path = self.path + (key, var_name) ds = self.blocks[block.name].add_dataset(dataset_struct, diff --git a/flopy/mf6/modflow/mfsimulation.py b/flopy/mf6/modflow/mfsimulation.py index 54827f9e1..86601a97a 100644 --- a/flopy/mf6/modflow/mfsimulation.py +++ b/flopy/mf6/modflow/mfsimulation.py @@ -353,6 +353,8 @@ class MFSimulation(PackageContainer): simulation name file _tdis_file simulation tdis file + sim_package_list + list of all "simulation level" packages Methods ------- @@ -599,7 +601,7 @@ def model_names(self): @classmethod def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe', sim_ws='.', strict=True, verbosity_level=1, load_only=None, - verify_data=True): + verify_data=False): """Load an existing model. Parameters @@ -801,8 +803,75 @@ def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe', ims_file.load(strict) instance.simulation_data.mfpath.set_last_accessed_path() + if verify_data: + instance.check() return instance + def check(self, f=None, verbose=True, level=1): + """ + Check model data for common errors. + + Parameters + ---------- + f : str or file handle + String defining file name or file handle for summary file + of check method output. If a string is passed a file handle + is created. If f is None, check method does not write + results to a summary file. (default is None) + verbose : bool + Boolean flag used to determine if check method results are + written to the screen + level : int + Check method analysis level. If level=0, summary checks are + performed. If level=1, full checks are performed. + + Returns + ------- + None + + Examples + -------- + + >>> import flopy + >>> m = flopy.modflow.Modflow.load('model.nam') + >>> m.check() + """ + # check instance for simulation-level check + chk_list = [] + + # check models + for model in self._models.values(): + print('Checking model "{}"...'.format(model.name)) + chk_list.append(model.check(f, verbose, level)) + + print('Checking for missing simulation packages...') + if self._tdis_file is None: + if chk_list: + chk_list[0]._add_to_summary( + 'Error', desc='\r No tdis package', package='model') + print('Error: no tdis package') + if len(self._ims_files) == 0: + if chk_list: + chk_list[0]._add_to_summary( + 'Error', desc='\r No solver package', package='model') + print('Error: no ims package') + return chk_list + + @property + def sim_package_list(self): + package_list = [] + if self._tdis_file is not None: + package_list.append(self._tdis_file) + for sim_package in self._ims_files.values(): + package_list.append(sim_package) + for sim_package in self._exchange_files.values(): + package_list.append(sim_package) + for sim_package in self._mover_files.values(): + package_list.append(sim_package) + for sim_package in self._other_files.values(): + package_list.append(sim_package) + return package_list + def load_package(self, ftype, fname, pname, strict, ref_path, dict_package_name=None, parent_package=None): """Load a package from a file. diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py index 77d7ccb73..044017189 100644 --- a/flopy/modflow/mf.py +++ b/flopy/modflow/mf.py @@ -244,7 +244,8 @@ def modeltime(self): 'tsmult': self.dis.tsmult.array} self._model_time = ModelTime(data_frame, self.dis.itmuni_dict[self.dis.itmuni], - self.dis.start_datetime, self.dis.steady) + self.dis.start_datetime, + self.dis.steady.array) return self._model_time @property @@ -278,7 +279,8 @@ def modelgrid(self): epsg=self._modelgrid.epsg, xoff=self._modelgrid.xoffset, yoff=self._modelgrid.yoffset, - angrot=self._modelgrid.angrot) + angrot=self._modelgrid.angrot, + nlay=self.dis.nlay) # resolve offsets xoff = self._modelgrid.xoffset diff --git a/flopy/modflow/mfbas.py b/flopy/modflow/mfbas.py index 74a06bc6e..8c0ecc766 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, check, get_neighbors +from ..utils import Util3d, get_neighbors class ModflowBas(Package): @@ -151,7 +151,7 @@ def __setattr__(self, key, value): else: super(ModflowBas, self).__setattr__(key, value) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check package data for common errors. @@ -181,7 +181,7 @@ def check(self, f=None, verbose=True, level=1): >>> m.bas6.check() """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) neighbors = get_neighbors(self.ibound.array) neighbors[ diff --git a/flopy/modflow/mfdis.py b/flopy/modflow/mfdis.py index 8d6393dc6..dbdc3d2f5 100644 --- a/flopy/modflow/mfdis.py +++ b/flopy/modflow/mfdis.py @@ -14,7 +14,7 @@ import numpy as np from ..pakbase import Package -from ..utils import Util2d, Util3d, check +from ..utils import Util2d, Util3d from ..utils.reference import SpatialReference, TemporalReference from ..utils.flopy_io import line_parse @@ -640,7 +640,7 @@ def write_file(self, check=True): f_dis.write(' {0:3s}\n'.format('TR')) f_dis.close() - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check dis package data for zero and negative thicknesses. @@ -669,7 +669,7 @@ def check(self, f=None, verbose=True, level=1): >>> m = flopy.modflow.Modflow.load('model.nam') >>> m.dis.check() """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) # make ibound of same shape as thicknesses/botm for quasi-3D models active = chk.get_active(include_cbd=True) diff --git a/flopy/modflow/mfmnw2.py b/flopy/modflow/mfmnw2.py index 5d55a34c6..83fd8e312 100644 --- a/flopy/modflow/mfmnw2.py +++ b/flopy/modflow/mfmnw2.py @@ -614,7 +614,7 @@ def sort_node_data(node_data): node_data = np.sort(node_data, order=['ztop'])[::-1] return node_data - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check mnw object for common errors. @@ -637,7 +637,7 @@ def check(self, f=None, verbose=True, level=1): chk : flopy.utils.check object """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) if self.losstype.lower() not in ['none', 'thiem', 'skin', 'general', 'sepecifycwc']: chk._add_to_summary(type='Error', k=self.k, i=self.i, j=self.j, @@ -646,6 +646,12 @@ def check(self, f=None, verbose=True, level=1): chk.summarize() return chk + def _get_check(self, f, verbose, level, checktype): + if checktype is not None: + return checktype(self, f=f, verbose=verbose, level=level) + else: + return check(self, f=f, verbose=verbose, level=level) + def _set_attributes_from_node_data(self): """ Populates the Mnw object attributes with values from node_data table. @@ -1300,7 +1306,7 @@ def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None): stress_period_data=stress_period_data, itmp=itmp, unitnumber=unitnumber, filenames=filenames) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check mnw2 package data for common errors. @@ -1329,7 +1335,7 @@ def check(self, f=None, verbose=True, level=1): >>> m = flopy.modflow.Modflow.load('model.nam') >>> m.mnw2.check() """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) # itmp if self.itmp[0] < 0: diff --git a/flopy/modflow/mfmnwi.py b/flopy/modflow/mfmnwi.py index 1aa5e47ca..84ea1e9df 100644 --- a/flopy/modflow/mfmnwi.py +++ b/flopy/modflow/mfmnwi.py @@ -1,9 +1,7 @@ import sys from ..utils.flopy_io import line_parse, pop_item - from ..pakbase import Package -from ..utils import check class ModflowMnwi(Package): @@ -247,7 +245,7 @@ def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None): extension='mnwi', unitnumber=unitnumber, filenames=filenames) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check mnwi package data for common errors. @@ -276,7 +274,7 @@ def check(self, f=None, verbose=True, level=1): >>> m = flopy.modflow.Modflow.load('model.nam') >>> m.mnwi.check() """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) if "MNW2" not in self.parent.get_package_list(): desc = '\r MNWI package present without MNW2 package.' chk._add_to_summary(type='Warning', value=0, diff --git a/flopy/modflow/mfoc.py b/flopy/modflow/mfoc.py index 57512238f..9d7e757d0 100644 --- a/flopy/modflow/mfoc.py +++ b/flopy/modflow/mfoc.py @@ -11,7 +11,6 @@ import sys from ..pakbase import Package -from ..utils import check class ModflowOc(Package): @@ -313,7 +312,7 @@ def __init__(self, model, \ self.parent.add_package(self) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check package data for common errors. @@ -343,7 +342,7 @@ def check(self, f=None, verbose=True, level=1): >>> m.oc.check() """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) dis = self.parent.get_package('DIS') if dis is None: dis = self.parent.get_package('DISU') diff --git a/flopy/modflow/mfrch.py b/flopy/modflow/mfrch.py index d82e3fe6a..1d4f09b47 100644 --- a/flopy/modflow/mfrch.py +++ b/flopy/modflow/mfrch.py @@ -11,7 +11,7 @@ import sys import numpy as np from ..pakbase import Package -from ..utils import Util2d, Transient2d, check +from ..utils import Util2d, Transient2d from ..modflow.mfparbc import ModflowParBc as mfparbc from ..utils.flopy_io import line_parse @@ -148,7 +148,8 @@ def __init__(self, model, nrchop=3, ipakcb=None, rech=1e-3, irch=0, self.np = 0 self.parent.add_package(self) - def check(self, f=None, verbose=True, level=1, RTmin=2e-8, RTmax=2e-4): + def check(self, f=None, verbose=True, level=1, RTmin=2e-8, RTmax=2e-4, + checktype=None): """ Check package data for common errors. @@ -182,7 +183,7 @@ def check(self, f=None, verbose=True, level=1, RTmin=2e-8, RTmax=2e-4): >>> m.rch.check() """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) if self.parent.bas6 is not None: active = self.parent.bas6.ibound.array.sum(axis=0) != 0 else: diff --git a/flopy/modflow/mfriv.py b/flopy/modflow/mfriv.py index 090f093ba..b27293da0 100644 --- a/flopy/modflow/mfriv.py +++ b/flopy/modflow/mfriv.py @@ -10,7 +10,7 @@ import sys import numpy as np from ..pakbase import Package -from ..utils import MfList, check +from ..utils import MfList from ..utils.recarray_utils import create_empty_recarray @@ -173,7 +173,7 @@ def __init__(self, model, ipakcb=None, stress_period_data=None, dtype=None, self.stress_period_data = MfList(self, stress_period_data) self.parent.add_package(self) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check package data for common errors. @@ -203,8 +203,9 @@ def check(self, f=None, verbose=True, level=1): >>> m.riv.check() """ - basechk = super(ModflowRiv, self).check(verbose=False) - chk = check(self, f=f, verbose=verbose, level=level) + basechk = super(ModflowRiv, self).check(verbose=False, + checktype=checktype) + chk = self._get_check(f, verbose, level, checktype) chk.summary_array = basechk.summary_array for per in self.stress_period_data.data.keys(): diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py index 85d452854..69d9cf6da 100644 --- a/flopy/modflow/mfsfr2.py +++ b/flopy/modflow/mfsfr2.py @@ -949,7 +949,7 @@ def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None): unit_number=unitnumber, filenames=filenames, options=options) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check sfr2 package data for common errors. diff --git a/flopy/modpath/mpsim.py b/flopy/modpath/mpsim.py index 45301b351..16917a122 100644 --- a/flopy/modpath/mpsim.py +++ b/flopy/modpath/mpsim.py @@ -9,7 +9,7 @@ """ import numpy as np from ..pakbase import Package -from ..utils import Util2d, Util3d, check +from ..utils import Util3d class ModpathSim(Package): @@ -128,7 +128,7 @@ def __init__(self, model, mp_name_file='mp.nam', mp_list_file='mp.list', self.parent.add_package(self) - def check(self, f=None, verbose=True, level=1): + def check(self, f=None, verbose=True, level=1, checktype=None): """ Check package data for common errors. @@ -153,7 +153,7 @@ def check(self, f=None, verbose=True, level=1): Examples -------- """ - chk = check(self, f=f, verbose=verbose, level=level) + chk = self._get_check(f, verbose, level, checktype) # MODPATH apparently produces no output if stoptime > last timepoint if self.options_dict['StopOption'] == 3 and self.options_dict[ diff --git a/flopy/mt3d/mt.py b/flopy/mt3d/mt.py index 0247759b6..73514a54c 100644 --- a/flopy/mt3d/mt.py +++ b/flopy/mt3d/mt.py @@ -341,7 +341,8 @@ def modeltime(self): self._model_time = ModelTime(data_frame, self.mf.dis.itmuni_dict[ self.mf.dis.itmuni], - self.dis.start_datetime, self.dis.steady) + self.dis.start_datetime, + self.dis.steady.array) return self._model_time @property @@ -355,11 +356,13 @@ def modelgrid(self): delr = self.btn.delr.array top = self.btn.htop.array botm = np.subtract(top, self.btn.dz.array.cumsum(axis=0)) + nlay = self.btn.nlay else: delc = self.mf.dis.delc.array delr = self.mf.dis.delr.array top = self.mf.dis.top.array botm = self.mf.dis.botm.array + nlay = self.mf.nlay if self.mf.bas6 is not None: ibound = self.mf.bas6.ibound.array else: @@ -374,7 +377,8 @@ def modelgrid(self): epsg=self._modelgrid.epsg, xoff=self._modelgrid.xoffset, yoff=self._modelgrid.yoffset, - angrot=self._modelgrid.angrot) + angrot=self._modelgrid.angrot, + nlay=nlay) # resolve offsets xoff = self._modelgrid.xoffset diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 95ef62b37..befc2b36a 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -25,35 +25,35 @@ class PackageInterface(object): @abc.abstractmethod def name(self): raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define name in child ' 'class to use this base class') @name.setter @abc.abstractmethod def name(self, name): raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define name in child ' 'class to use this base class') @property @abc.abstractmethod def parent(self): raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define parent in child ' 'class to use this base class') @parent.setter @abc.abstractmethod def parent(self, name): raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define parent in child ' 'class to use this base class') @property @abc.abstractmethod def package_type(self): raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define package_type in child ' 'class to use this base class') @property @@ -61,13 +61,13 @@ def package_type(self): def data_list(self): # [data_object, data_object, ...] raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define data_list in child ' 'class to use this base class') @abc.abstractmethod def export(self, f, **kwargs): raise NotImplementedError( - 'must define get_model_dim_arrays in child ' + 'must define export in child ' 'class to use this base class') @property @@ -77,6 +77,270 @@ def plotable(self): 'must define plotable in child ' 'class to use this base class') + @property + def has_stress_period_data(self): + return self.__dict__.get('stress_period_data', None) is not None + + @staticmethod + def _check_thresholds(chk, array, active, thresholds, name): + """Checks array against min and max threshold values.""" + mn, mx = thresholds + chk.values(array, active & (array < mn), + '{} values below checker threshold of {}' + .format(name, mn), 'Warning') + chk.values(array, active & (array > mx), + '{} values above checker threshold of {}' + .format(name, mx), 'Warning') + + @staticmethod + def _confined_layer_check(chk): + return + + def _other_xpf_checks(self, chk, active): + # check for negative hani + chk.values(self.__dict__['hani'].array, + active & (self.__dict__['hani'].array < 0), + 'negative horizontal anisotropy values', 'Error') + + # check vkcb if there are any quasi-3D layers + if self.parent.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: + # assign 1 instead of zero as default value that + # won't violate checker + # (allows for same structure as other checks) + vkcb[l, :, :] = 1 + chk.values(vkcb, active & (vkcb <= 0), + 'zero or negative quasi-3D confining bed Kv values', + 'Error') + self._check_thresholds(chk, vkcb, active, + chk.property_threshold_values['vkcb'], + 'quasi-3D confining bed Kv') + + @staticmethod + def _get_nan_exclusion_list(): + return [] + + def _get_check(self, f, verbose, level, checktype): + if checktype is not None: + return checktype(self, f=f, verbose=verbose, level=level) + else: + return check(self, f=f, verbose=verbose, level=level) + + def _check_oc(self, f=None, verbose=True, level=1, checktype=None): + spd_inds_valid = True + chk = self._get_check(f, verbose, level, checktype) + spd = getattr(self, 'stress_period_data') + nan_exclusion_list = self._get_nan_exclusion_list() + for per in spd.data.keys(): + if isinstance(spd.data[per], np.recarray): + spdata = self.stress_period_data.data[per] + inds = chk._get_cell_inds(spdata) + + # General BC checks + # check for valid cell indices + spd_inds_valid = \ + chk._stress_period_data_valid_indices(spdata) + + # first check for and list nan values + chk._stress_period_data_nans(spdata, nan_exclusion_list) + + if spd_inds_valid: + # next check for BCs in inactive cells + chk._stress_period_data_inactivecells(spdata) + + # More specific BC checks + # check elevations in the ghb, drain, and riv packages + if self.name[0] in check.bc_stage_names.keys(): + # check that bc elevations are above model + # cell bottoms -- also checks for nan values + elev_name = chk.bc_stage_names[self.name[0]] + mg = self.parent.modelgrid + botms = mg.botm[inds] + test = spdata[elev_name] < botms + en = 'BC elevation below cell bottom' + chk.stress_period_data_values(spdata, + test, + col=elev_name, + error_name=en, + error_type='Error') + + chk.summarize() + return chk + + def _get_kparams(self): + # build model specific parameter lists + kparams_all = {'hk': 'horizontal hydraulic conductivity', + 'vka': 'vertical hydraulic conductivity', + 'k': 'horizontal hydraulic conductivity', + 'k22': 'hydraulic conductivity second axis', + 'k33': 'vertical hydraulic conductivity'} + kparams = {} + vka_param = None + for kp, name in kparams_all.items(): + if kp in self.__dict__: + kparams[kp] = name + if 'hk' in self.__dict__: + 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 + vka_param = kparams.pop('vka') + elif 'k33' in self.__dict__: + vka = self.k33.array + vka_param = kparams.pop('k33') + else: + vka = None + if vka is not None: + vka = vka.copy() + return kparams, hk, vka, vka_param + + def _check_flowp(self, f=None, verbose=True, level=1, checktype=None): + chk = self._get_check(f, verbose, level, checktype) + active = chk.get_active() + + # build model specific parameter lists + kparams, hk, vka, vka_param = self._get_kparams() + + # check for zero or negative values of hydraulic conductivity, + # anisotropy, and quasi-3D confining beds + for kp, name in kparams.items(): + if self.__dict__[kp].array is not None: + chk.values(self.__dict__[kp].array, + active & (self.__dict__[kp].array <= 0), + 'zero or negative {} values'.format(name), + 'Error') + + if 'hani' in self.__dict__: + self._other_xpf_checks(chk, active) + + # check for unusually high or low values of hydraulic conductivity + # convert vertical anisotropy to Kv for checking + if vka is not None: + if 'layvka' in self.__dict__: + for l in range(vka.shape[0]): + vka[l] *= hk[l] if self.layvka.array[l] != 0 else 1 + self._check_thresholds(chk, vka, active, + chk.property_threshold_values['vka'], + vka_param) + + for kp, name in kparams.items(): + if self.__dict__[kp].array is not None: + self._check_thresholds(chk, self.__dict__[kp].array, + active, + chk.property_threshold_values[kp], + name) + if self.name[0] in ['UPW', 'LPF']: + storage_coeff = 'STORAGECOEFFICIENT' in self.options or \ + ('storagecoefficient' in self.__dict__ and + self.storagecoefficient.get_data()) + self._check_storage(chk, storage_coeff) + chk.summarize() + return chk + + def check(self, f=None, verbose=True, level=1, checktype=None): + """ + Check package data for common errors. + + Parameters + ---------- + f : str or file handle + String defining file name or file handle for summary file + of check method output. If a sting is passed a file handle + is created. If f is None, check method does not write + results to a summary file. (default is None) + verbose : bool + Boolean flag used to determine if check method results are + written to the screen + level : int + Check method analysis level. If level=0, summary checks are + performed. If level=1, full checks are performed. + checktype : check + Checker type to be used. By default class check is used from + check.py. + + Returns + ------- + None + + Examples + -------- + + >>> import flopy + >>> m = flopy.modflow.Modflow.load('model.nam') + >>> m.dis.check() + + """ + chk = None + + if self.has_stress_period_data and self.name[0] != 'OC' and \ + self.package_type.upper() != 'OC': + chk = self._check_oc(f, verbose, level, checktype) + # check property values in upw and lpf packages + elif self.name[0] in ['UPW', 'LPF'] or \ + self.package_type.upper() in ['NPF']: + chk = self._check_flowp(f, verbose, level, checktype) + elif self.package_type.upper() in ['STO']: + chk = self._get_check(f, verbose, level, checktype) + storage_coeff = self.storagecoefficient.get_data() + if storage_coeff is None: + storage_coeff = False + self._check_storage(chk, storage_coeff) + else: + txt = 'check method not implemented for ' + \ + '{} Package.'.format(self.name[0]) + if f is not None: + if isinstance(f, str): + pth = os.path.join(self.parent.model_ws, f) + f = open(pth, 'w') + f.write(txt) + f.close() + if verbose: + print(txt) + return chk + + def _check_storage(self, chk, storage_coeff): + # only check storage if model is transient + if not np.all(self.parent.modeltime.steady_state): + active = chk.get_active() + # do the same for storage if the model is transient + sarrays = {'ss': self.ss.array, 'sy': self.sy.array} + # convert to specific for checking + if storage_coeff: + desc = '\r STORAGECOEFFICIENT option is ' + \ + 'activated, storage values are read ' + \ + 'storage coefficients' + chk._add_to_summary(type='Warning', desc=desc) + + chk.values(sarrays['ss'], active & (sarrays['ss'] < 0), + 'zero or negative specific storage values', 'Error') + self._check_thresholds(chk, sarrays['ss'], active, + chk.property_threshold_values['ss'], + 'specific storage') + + # only check specific yield for convertible layers + if 'laytyp' in self.__dict__: + inds = np.array( + [True if l > 0 or l < 0 and 'THICKSTRT' in self.options + else False for l in self.laytyp]) + sarrays['sy'] = sarrays['sy'][inds, :, :] + active = active[inds, :, :] + else: + iconvert = self.iconvert.array + for ishape in np.ndindex(active.shape): + if active[ishape]: + active[ishape] = iconvert[ishape] > 0 or \ + iconvert[ishape] < 0 + chk.values(sarrays['sy'], active & (sarrays['sy'] < 0), + 'zero or negative specific yield values', 'Error') + self._check_thresholds(chk, sarrays['sy'], active, + chk.property_threshold_values['sy'], + 'specific yield') + class Package(PackageInterface): """ @@ -315,205 +579,21 @@ def get_sfac_columns(): """ return [] - def check(self, f=None, verbose=True, level=1): - """ - Check package data for common errors. - - Parameters - ---------- - f : str or file handle - String defining file name or file handle for summary file - of check method output. If a sting is passed a file handle - is created. If f is None, check method does not write - results to a summary file. (default is None) - verbose : bool - Boolean flag used to determine if check method results are - written to the screen - level : int - Check method analysis level. If level=0, summary checks are - performed. If level=1, full checks are performed. - - Returns - ------- - None - - Examples - -------- - - >>> import flopy - >>> m = flopy.modflow.Modflow.load('model.nam') - >>> m.dis.check() - - """ - chk = None - - if self.__dict__.get('stress_period_data', None) is not None \ - and self.name[0] != 'OC': - spd_inds_valid = True - chk = check(self, f=f, verbose=verbose, level=level) - spd = getattr(self, 'stress_period_data') - for per in spd.data.keys(): - if isinstance(spd.data[per], np.recarray): - spdata = self.stress_period_data.data[per] - inds = (spdata.k, spdata.i, spdata.j) \ - if self.parent.structured \ - else (spdata.node) - - # General BC checks - # check for valid cell indices - spd_inds_valid = \ - chk._stress_period_data_valid_indices(spdata) - - # first check for and list nan values - chk._stress_period_data_nans(spdata) - - if spd_inds_valid: - # next check for BCs in inactive cells - chk._stress_period_data_inactivecells(spdata) - - # More specific BC checks - # check elevations in the ghb, drain, and riv packages - if self.name[0] in check.bc_stage_names.keys(): - # check that bc elevations are above model - # cell bottoms -- also checks for nan values - elev_name = chk.bc_stage_names[self.name[0]] - botms = self.parent.dis.botm.array[inds] - test = spdata[elev_name] < botms - en = 'BC elevation below cell bottom' - chk.stress_period_data_values(spdata, - test, - col=elev_name, - error_name=en, - error_type='Error') - - chk.summarize() - - # check property values in upw and lpf packages - elif self.name[0] in ['UPW', 'LPF']: - - chk = check(self, f=f, verbose=verbose, level=level) - active = chk.get_active() - - # check for confined layers above convertible layers - confined = False - thickstrt = False - for option in self.options: - if option.lower() == 'thickstrt': - thickstrt = True - for i, l in enumerate(self.laytyp.array.tolist()): - if l == 0 or l < 0 and thickstrt: - confined = True - continue - if confined and l > 0: - desc = '\r LAYTYP: unconfined (convertible) ' + \ - 'layer below confined layer' - chk._add_to_summary(type='Warning', desc=desc) - - # check for zero or negative values of hydraulic conductivity, - # anisotropy, and quasi-3D confining beds - kparams = {'hk': 'horizontal hydraulic conductivity', - 'vka': 'vertical hydraulic conductivity'} - for kp, name in kparams.items(): - chk.values(self.__dict__[kp].array, - active & (self.__dict__[kp].array <= 0), - 'zero or negative {} values'.format(name), 'Error') - - # check for negative hani - chk.values(self.__dict__['hani'].array, - active & (self.__dict__['hani'].array < 0), - 'negative horizontal anisotropy values', 'Error') - - def check_thresholds(array, active, thresholds, name): - """Checks array against min and max threshold values.""" - mn, mx = thresholds - chk.values(array, active & (array < mn), - '{} values below checker threshold of {}' - .format(name, mn), 'Warning') - chk.values(array, active & (array > mx), - '{} values above checker threshold of {}' - .format(name, mx), 'Warning') - - # check for unusually high or low values of hydraulic conductivity - # convert vertical anisotropy to Kv for checking - if self.layvka.sum() > 0: - vka = self.vka.array.copy() - for l in range(vka.shape[0]): - vka[l] *= self.hk.array[l] if self.layvka.array[ - l] != 0 else 1 - check_thresholds(vka, active, - chk.property_threshold_values['vka'], - kparams.pop('vka')) - - for kp, name in kparams.items(): - check_thresholds(self.__dict__[kp].array, active, - chk.property_threshold_values[kp], - name) - - # check vkcb if there are any quasi-3D layers - if self.parent.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: - # assign 1 instead of zero as default value that - # won't violate checker - # (allows for same structure as other checks) - vkcb[l, :, :] = 1 - chk.values(vkcb, active & (vkcb <= 0), - 'zero or negative quasi-3D confining bed Kv values', - 'Error') - check_thresholds(vkcb, active, - chk.property_threshold_values['vkcb'], - 'quasi-3D confining bed Kv') - - # only check storage if model is transient - if not np.all(self.parent.dis.steady.array): - - # do the same for storage if the model is transient - sarrays = {'ss': self.ss.array, 'sy': self.sy.array} - # convert to specific for checking - if 'STORAGECOEFFICIENT' in self.options: - desc = '\r STORAGECOEFFICIENT option is ' + \ - 'activated, storage values are read ' + \ - 'storage coefficients' - chk._add_to_summary(type='Warning', desc=desc) - tshape = (self.parent.nlay, self.parent.nrow, - self.parent.ncol) - sarrays['ss'].shape != tshape - sarrays['sy'].shape != tshape - - chk.values(sarrays['ss'], active & (sarrays['ss'] < 0), - 'zero or negative specific storage values', 'Error') - check_thresholds(sarrays['ss'], active, - chk.property_threshold_values['ss'], - 'specific storage') - - # only check specific yield for convertible layers - inds = np.array( - [True if l > 0 or l < 0 and 'THICKSTRT' in self.options - else False for l in self.laytyp]) - sarrays['sy'] = sarrays['sy'][inds, :, :] - active = active[inds, :, :] - chk.values(sarrays['sy'], active & (sarrays['sy'] < 0), - 'zero or negative specific yield values', 'Error') - check_thresholds(sarrays['sy'], active, - chk.property_threshold_values['sy'], - 'specific yield') - chk.summarize() - - else: - txt = 'check method not implemented for ' + \ - '{} Package.'.format(self.name[0]) - if f is not None: - if isinstance(f, str): - pth = os.path.join(self.parent.model_ws, f) - f = open(pth, 'w') - f.write(txt) - f.close() - if verbose: - print(txt) - return chk + def _confined_layer_check(self, chk): + # check for confined layers above convertible layers + confined = False + thickstrt = False + for option in self.options: + if option.lower() == 'thickstrt': + thickstrt = True + for i, l in enumerate(self.laytyp.array.tolist()): + if l == 0 or l < 0 and thickstrt: + confined = True + continue + if confined and l > 0: + desc = '\r LAYTYP: unconfined (convertible) ' + \ + 'layer below confined layer' + chk._add_to_summary(type='Warning', desc=desc) def level1_arraylist(self, idx, v, name, txt): ndim = v.ndim diff --git a/flopy/seawat/swt.py b/flopy/seawat/swt.py index d2b88dece..410b63652 100644 --- a/flopy/seawat/swt.py +++ b/flopy/seawat/swt.py @@ -155,7 +155,8 @@ def modeltime(self): 'tsmult': self.dis.tsmult.array} self._model_time = ModelTime(data_frame, self.dis.itmuni_dict[self.dis.itmuni], - self.dis.start_datetime, self.dis.steady) + self.dis.start_datetime, + self.dis.steady.array) return self._model_time @property @@ -178,7 +179,8 @@ def modelgrid(self): epsg=self._modelgrid.epsg, xoff=self._modelgrid.xoffset, yoff=self._modelgrid.yoffset, - angrot=self._modelgrid.angrot) + angrot=self._modelgrid.angrot, + nlay=self.dis.nlay) # resolve offsets xoff = self._modelgrid.xoffset diff --git a/flopy/utils/check.py b/flopy/utils/check.py index 550150680..d14aeda00 100644 --- a/flopy/utils/check.py +++ b/flopy/utils/check.py @@ -61,9 +61,12 @@ class check: package_check_levels = {'sfr': 1} property_threshold_values = {'hk': (1e-11, 1e5), + 'k': (1e-11, 1e5), + 'k22': (1e-11, 1e5), # after Schwartz and Zhang, table 4.4 'hani': None, 'vka': (1e-11, 1e5), + 'k33': (1e-11, 1e5), 'vkcb': (1e-11, 1e5), 'ss': (1e-6, 1e-2), 'sy': (0.01, 0.5)} @@ -90,7 +93,10 @@ def __init__(self, package, f=None, verbose=True, level=1, self.prefix = '{} MODEL DATA VALIDATION SUMMARY'.format( self.model.name) self.package = package - self.structured = self.model.structured + if 'structured' in self.model.__dict__: + self.structured = self.model.structured + else: + self.structured = (self.model.modelgrid.grid_type == 'structured') self.verbose = verbose self.level = level self.passed = [] @@ -207,24 +213,7 @@ def _boolean_compare(self, array, col1, col2, return txt def _get_summary_array(self, array=None): - - if self.structured: - # include node column for structured grids (useful for indexing) - dtype = np.dtype([('type', np.object), - ('package', np.object), - ('k', np.int), - ('i', np.int), - ('j', np.int), - ('value', np.float), - ('desc', np.object) - ]) - else: - dtype = np.dtype([('type', np.object), - ('package', np.object), - ('node', np.int), - ('value', np.float), - ('desc', np.object) - ]) + dtype = self._get_dtype() if array is None: return np.recarray((0), dtype=dtype) ra = recarray(array, dtype) @@ -249,25 +238,10 @@ def _txt_footer(self, headertxt, txt, testname, passed=False, def _stress_period_data_valid_indices(self, stress_period_data): """Check that stress period data inds are valid for model grid.""" - spd_inds_valid = True - if self.model.has_package('DIS') and \ - {'k', 'i', 'j'}.intersection( - set(stress_period_data.dtype.names)) != {'k', 'i', 'j'}: - self._add_to_summary(type='Error', - desc='\r Stress period data missing k, i, j for structured grid.') - spd_inds_valid = False - elif self.model.has_package('DISU') and \ - 'node' not in stress_period_data.dtype.names: - self._add_to_summary(type='Error', - desc='\r Stress period data missing node number for unstructured grid.') - spd_inds_valid = False + spd_inds_valid = self._has_cell_indices(stress_period_data) # check for BCs indices that are invalid for grid - inds = (stress_period_data.k, - stress_period_data.i, - stress_period_data.j) if self.structured else ( - stress_period_data.node) - + inds = self._get_cell_inds(stress_period_data) isvalid = self.isvalid(inds) if not np.all(isvalid): sa = self._list_spd_check_violations(stress_period_data, ~isvalid, @@ -281,12 +255,13 @@ def _stress_period_data_valid_indices(self, stress_period_data): self.append_passed('BC indices valid') return spd_inds_valid - def _stress_period_data_nans(self, stress_period_data): + def _stress_period_data_nans(self, stress_period_data, nan_excl_list): """Check for and list any nans in stress period data.""" isnan = np.array([np.isnan(stress_period_data[c]) for c in stress_period_data.dtype.names - if not isinstance(stress_period_data.dtype[c], - np.object)]).transpose() + if not (stress_period_data.dtype[c].name + == 'object') and c not in + nan_excl_list]).transpose() if np.any(isnan): row_has_nan = np.any(isnan, axis=1) sa = self._list_spd_check_violations(stress_period_data, @@ -302,10 +277,12 @@ def _stress_period_data_nans(self, stress_period_data): def _stress_period_data_inactivecells(self, stress_period_data): """Check for and list any stress period data in cells with ibound=0.""" spd = stress_period_data - inds = (spd.k, spd.i, spd.j) if self.structured else (spd.node) + inds = self._get_cell_inds(spd) msg = 'BC in inactive cell' - if 'BAS6' in self.model.get_package_list(): - ibnd = self.package.parent.bas6.ibound.array[inds] + + idomain = self.model.modelgrid.idomain + if idomain is not None: + ibnd = idomain[inds] if np.any(ibnd == 0): sa = self._list_spd_check_violations(stress_period_data, @@ -326,15 +303,12 @@ def _list_spd_check_violations(self, stress_period_data, criteria, name, k,i,j indices, values, and description of error for each row in stress_period_data where criteria=True. """ - inds_col = ['k', 'i', 'j'] if self.structured else ['node'] + inds_col = self._get_cell_inds_names() # inds = stress_period_data[criteria][inds_col]\ # .reshape(stress_period_data[criteria].shape + (-1,)) # inds = np.atleast_2d(np.squeeze(inds.tolist())) inds = stress_period_data[criteria] - a = inds[inds_col[0]] - if len(inds_col) > 1: - for n in inds_col[1:]: - a = np.concatenate((a, inds[n])) + a = self._get_cellid_cols(inds, inds_col) inds = a.view(int) inds = inds.reshape(stress_period_data[criteria].shape + (-1,)) @@ -347,6 +321,14 @@ def _list_spd_check_violations(self, stress_period_data, criteria, tp = [error_type] * len(v) return self._get_summary_array(np.column_stack([tp, pn, inds, v, en])) + @staticmethod + def _get_cellid_cols(inds, inds_col): + a = inds[inds_col[0]] + if len(inds_col) > 1: + for n in inds_col[1:]: + a = np.concatenate((a, inds[n])) + return a + def append_passed(self, message): """Add a check to the passed list if it isn't already in there.""" self.passed.append(message) if message not in self.passed else None @@ -371,16 +353,18 @@ def isvalid(self, inds): if isinstance(inds, np.ndarray): inds = [inds] - if 'DIS' in self.model.get_package_list() and len(inds) == 3: - dis = self.model.dis - k = inds[0] < dis.nlay - i = inds[1] < dis.nrow - j = inds[2] < dis.ncol - return k | i | j - - elif 'DISU' in self.model.get_package_list() and len(inds) == 1: - return inds < self.model.disu.nodes - + mg = self.model.modelgrid + if mg.grid_type == 'structured' and len(inds) == 3: + k = inds[0] < mg.nlay + i = inds[1] < mg.nrow + j = inds[2] < mg.ncol + return k & i & j + elif mg.grid_type == 'vertex' and len(inds) == 2: + lay = inds[0] < mg.nlay + cpl = inds[1] < mg.ncpl + return lay & cpl + elif mg.grid_type == 'unstructured' and len(inds) == 1: + return inds[0] < mg.nnodes else: return np.zeros(inds[0].shape, dtype=bool) @@ -399,15 +383,21 @@ def get_active(self, include_cbd=False): active : 3-D boolean array True where active. """ - if 'DIS' in self.model.get_package_list(): - dis = self.model.dis - inds = (dis.nlay, dis.nrow, dis.ncol) + mg = self.model.modelgrid + if mg.grid_type == 'structured': + inds = (mg.nlay_nocbd, mg.nrow, mg.ncol) + elif mg.grid_type == 'vertex': + inds = (mg.nlay, mg.ncpl) else: - dis = self.model.disu - inds = dis.nodes + inds = mg.nnodes include_cbd = False if 'BAS6' in self.model.get_package_list(): + if 'DIS' in self.model.get_package_list(): + dis = self.model.dis + else: + dis = self.model.disu + # make ibound of same shape as thicknesses/botm for quasi-3D models if include_cbd and dis.laycbd.sum() > 0: ncbd = np.sum(dis.laycbd.array > 0) @@ -557,6 +547,48 @@ def summarize(self): if self.f is not None: print(' see {} for details.\n'.format(self.summaryfile)) + # start of older model specific code + def _has_cell_indices(self, stress_period_data): + if self.model.has_package('DIS') and \ + {'k', 'i', 'j'}.intersection( + set(stress_period_data.dtype.names)) != {'k', 'i', 'j'}: + self._add_to_summary(type='Error', + desc='\r Stress period data missing k, ' + 'i, j for structured grid.') + return False + elif self.model.has_package('DISU') and \ + 'node' not in stress_period_data.dtype.names: + self._add_to_summary(type='Error', + desc='\r Stress period data missing ' + 'node number for unstructured grid.') + return False + return True + + def _get_cell_inds(self, spd): + return (spd.k, spd.i, spd.j) if self.structured else (spd.node) + + def _get_cell_inds_names(self): + return ['k', 'i', 'j'] if self.structured else ['node'] + + def _get_dtype(self): + if self.structured: + # include node column for structured grids (useful for indexing) + return np.dtype([('type', np.object), + ('package', np.object), + ('k', np.int), + ('i', np.int), + ('j', np.int), + ('value', np.float), + ('desc', np.object) + ]) + else: + return np.dtype([('type', np.object), + ('package', np.object), + ('node', np.int), + ('value', np.float), + ('desc', np.object) + ]) + def _fmt_string_list(array, float_format='{}'): fmt_string = [] @@ -616,7 +648,8 @@ def _print_rec_array(array, cols=None, delimiter=' ', float_format='{:.6f}'): def fields_view(arr, fields): """ creates view of array that only contains the fields in fields. - http://stackoverflow.com/questions/15182381/how-to-return-a-view-of-several-columns-in-numpy-structured-array + http://stackoverflow.com/questions/15182381/how-to-return-a-view-of- + several-columns-in-numpy-structured-array """ dtype2 = np.dtype({name: arr.dtype.fields[name] for name in fields}) return np.ndarray(arr.shape, dtype2, arr, 0, arr.strides) @@ -635,7 +668,8 @@ def get_neighbors(a): ------- 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. + 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 @@ -649,3 +683,95 @@ def get_neighbors(a): 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, package, f=None, verbose=True, level=1, + property_threshold_values={}): + super(mf6check, self).__init__(package, f, verbose, level, + property_threshold_values) + if hasattr(package, 'model_or_sim'): + self.model = package.model_or_sim + + @staticmethod + def _get_cellid_cols(inds, inds_col): + a = inds[inds_col[0]] + return np.asarray(a.tolist()) + + + def _get_cell_inds(self, spd): + hnames = () + if 'cellid' in spd.dtype.names: + cellid = spd.cellid + elif 'cellid1' in spd.dtype.names: + cellid = spd.cellid1 + else: + return None + + for item in zip(*cellid): + hnames += (np.ndarray(shape=(len(item),), + buffer=np.array(item), dtype=np.int32),) + return hnames + + def _get_dtype(self): + mg = self.model.modelgrid + if mg.grid_type == 'structured': + return np.dtype([('type', np.object), + ('package', np.object), + ('k', np.int), + ('i', np.int), + ('j', np.int), + ('value', np.float), + ('desc', np.object) + ]) + elif mg.grid_type == 'vertex': + return np.dtype([('type', np.object), + ('package', np.object), + ('lay', np.int), + ('cell', np.int), + ('value', np.float), + ('desc', np.object) + ]) + else: + return np.dtype([('type', np.object), + ('package', np.object), + ('node', np.int), + ('value', np.float), + ('desc', np.object) + ]) + + def _has_cell_indices(self, stress_period_data): + mg = self.model.modelgrid + if mg.grid_type == 'structured' or mg.grid_type == 'vertex' or \ + mg.grid_type == 'unstructured': + if 'cellid' not in set(stress_period_data.dtype.names) and \ + 'cellid1' not in set(stress_period_data.dtype.names): + self._add_to_summary(type='Error', + desc='\r Stress period data missing ' + 'cellid.') + return False + return True + + def _get_cell_inds_names(self): + return ['cellid'] + + def get_active(self, include_cbd=False): + """Returns a boolean array of active cells for the model. + + Parameters + ---------- + include_cbd : boolean + Does not apply to MF6 models, always false. + + Returns + ------- + active : 3-D boolean array + True where active. + """ + mg = self.model.modelgrid + idomain = mg.idomain + if idomain is None: + return np.ones(shape=mg.shape, dtype=bool) + else: + id_active_zero = idomain - np.ones(shape=mg.shape, dtype=np.int32) + return np.invert(np.ma.make_mask(id_active_zero, shrink=False))