From 6d29e387ea5085c9b532d0ee6586f346355cf089 Mon Sep 17 00:00:00 2001 From: Scott Paulinski Date: Fri, 5 Nov 2021 10:54:27 -0700 Subject: [PATCH 1/6] fix(array): Getting array data (.array) fixed so that empty stress periods repeat data from the last previous stress period with data. If there is no previous stress period with data, None is used for MFDataList arrays and an array of zeros is used for MFDataArray arrays. --- flopy/mf6/data/mfdataarray.py | 71 ++++++++++++++++++++++------------- flopy/mf6/data/mfdatalist.py | 8 ++-- 2 files changed, 48 insertions(+), 31 deletions(-) diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index 1c7e58a9c4..318c46d991 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -1599,6 +1599,49 @@ def store_internal( check_data, ) + def _get_array(self, num_sp, apply_mult, **kwargs): + """Returns all data up to stress period num_sp in a single array. + If `layer` is None, returns all data for time `layer`. + + Parameters + ---------- + num_sp : int + Zero-based stress period of data to return all data up to + apply_mult : bool + Whether to apply multiplier to data prior to returning it + + """ + data = None + output = None + for sp in range(0, num_sp): + if sp in self._data_storage: + self.get_data_prep(sp) + data = super().get_data(apply_mult=apply_mult, **kwargs) + data = np.expand_dims(data, 0) + else: + # if there is no previous data provide array of + # zeros, otherwise provide last array of data found + if data is None: + # get any data + data_dimensions = None + for key in self._data_storage.keys(): + self.get_data_prep(key) + data = super().get_data( + apply_mult=apply_mult, **kwargs + ) + break + if data is not None: + if self.structure.type == DatumType.integer: + data = np.full_like(data, 0) + else: + data = np.full_like(data, 0.0) + data = np.expand_dims(data, 0) + if output is None or data is None: + output = data + else: + output = np.concatenate((output, data)) + return output + def get_data(self, layer=None, apply_mult=True, **kwargs): """Returns the data associated with stress period key `layer`. If `layer` is None, returns all data for time `layer`. @@ -1613,38 +1656,14 @@ def get_data(self, layer=None, apply_mult=True, **kwargs): """ if self._data_storage is not None and len(self._data_storage) > 0: if layer is None: - output = None sim_time = self._data_dimensions.package_dim.model_dim[ 0 ].simulation_time num_sp = sim_time.get_num_stress_periods() if "array" in kwargs: - data = None - for sp in range(0, num_sp): - if sp in self._data_storage: - self.get_data_prep(sp) - data = super().get_data( - apply_mult=apply_mult, **kwargs - ) - data = np.expand_dims(data, 0) - else: - if data is None: - # get any data - self.get_data_prep(self._data_storage.key()[0]) - data = super().get_data( - apply_mult=apply_mult, **kwargs - ) - data = np.expand_dims(data, 0) - if self.structure.type == DatumType.integer: - data = np.full_like(data, 0) - else: - data = np.full_like(data, 0.0) - if output is None: - output = data - else: - output = np.concatenate((output, data)) - return output + return self._get_array(num_sp, apply_mult, **kwargs) else: + output = None for sp in range(0, num_sp): data = None if sp in self._data_storage: diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index 34308cbcd8..400e157408 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -1641,14 +1641,12 @@ def get_data(self, key=None, apply_mult=False, **kwargs): 0 ].simulation_time num_sp = sim_time.get_num_stress_periods() + data = None for sp in range(0, num_sp): if sp in self._data_storage: self.get_data_prep(sp) - output.append( - super().get_data(apply_mult=apply_mult) - ) - else: - output.append(None) + data = super().get_data(apply_mult=apply_mult) + output.append(data) return output else: output = {} From 62b1b5a3b2c27d0ae24c28d5c32cf78b7490c62e Mon Sep 17 00:00:00 2001 From: Scott Paulinski Date: Fri, 5 Nov 2021 10:56:18 -0700 Subject: [PATCH 2/6] fix(array): added test cases for getting array data (.array) --- autotest/t505_test.py | 62 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/autotest/t505_test.py b/autotest/t505_test.py index ac53ac37aa..81193ed875 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -352,6 +352,67 @@ def test_multi_model(): sim.run_simulation() +def test_array(): + import flopy + + mf6 = flopy.mf6 + sim_name = "testsim" + model_name = "testmodel" + out_dir = "tmp" + + tdis_name = "{}.tdis".format(sim_name) + sim = mf6.MFSimulation( + sim_name=sim_name, version="mf6", exe_name=exe_name, sim_ws=out_dir + ) + tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0), (6.0, 3, 1.0), (6.0, 3, 1.0)] + tdis = mf6.ModflowTdis(sim, time_units="DAYS", nper=4, perioddata=tdis_rc) + + model = mf6.ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + + dis = mf6.ModflowGwfdis( + model, + length_units="FEET", + nlay=1, + nrow=2, + ncol=2, + delr=500.0, + delc=500.0, + top=100.0, + botm=50.0, + filename="{}.dis".format(model_name), + ) + irch = {1: [[0, 1], [2, 1]], 2: [[0, 1], [2, 3]]} + rcha = mf6.ModflowGwfrcha(model, irch=irch, recharge={0: 1.0, 2: 2.0}) + val_irch = rcha.irch.array.sum(axis=(1, 2, 3)) + assert val_irch[0] == 0 + assert val_irch[1] == 4 + assert val_irch[2] == 6 + assert val_irch[3] == 6 + val_rch = rcha.recharge.array.sum(axis=(1, 2, 3)) + assert val_rch[0] == 4.0 + assert val_rch[1] == 4.0 + assert val_rch[2] == 8.0 + assert val_rch[3] == 8.0 + + welspdict = {1: [[(0, 0, 0), -25.0, 0.0]], 2: [[(0, 0, 0), 25.0, 0.0]]} + wel = flopy.mf6.ModflowGwfwel( + model, + print_input=True, + print_flows=True, + stress_period_data=welspdict, + save_flows=False, + auxiliary="CONCENTRATION", + pname="WEL-1", + ) + wel_array = wel.stress_period_data.array + assert wel_array[0] is None + assert wel_array[1][0][1] == -25.0 + assert wel_array[2][0][1] == 25.0 + assert wel_array[3][0][1] == 25.0 + + def test_np001(): # init paths test_ex_name = "np001" @@ -3580,6 +3641,7 @@ def test_transport(): if __name__ == "__main__": + test_array() test_multi_model() test_np001() test_np002() From 4c486837848115f5e38964c145c28e5bccc115b7 Mon Sep 17 00:00:00 2001 From: Scott Paulinski Date: Mon, 29 Nov 2021 12:49:09 -0800 Subject: [PATCH 3/6] fix(array and getdata): array and getdata now return appropriate values for "nodata" and "repeating last stress period data" cases --- flopy/mf6/data/mfdata.py | 10 +++ flopy/mf6/data/mfdataarray.py | 29 ++++++- flopy/mf6/data/mfdatalist.py | 28 ++++++- flopy/mf6/data/mfdatascalar.py | 13 ++- flopy/mf6/data/mfdatastorage.py | 15 +++- flopy/mf6/mfpackage.py | 144 +++++++++++++++++++++----------- 6 files changed, 175 insertions(+), 64 deletions(-) diff --git a/flopy/mf6/data/mfdata.py b/flopy/mf6/data/mfdata.py index 4a49eafa39..9f5e2a0193 100644 --- a/flopy/mf6/data/mfdata.py +++ b/flopy/mf6/data/mfdata.py @@ -100,6 +100,12 @@ def get_data_prep(self, transient_key=0): self._verify_sp(transient_key) self._current_key = transient_key if transient_key not in self._data_storage: + if ( + isinstance(transient_key, tuple) + and transient_key[0] in self._data_storage + ): + self._current_key = transient_key[0] + return self.add_transient_key(transient_key) def _set_data_prep(self, data, transient_key=0): @@ -268,6 +274,10 @@ def __repr__(self): def __str__(self): return str(self._get_storage_obj()) + @property + def path(self): + return self._path + @property def array(self): kwargs = {"array": True} diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index 318c46d991..0ff3fbc9f8 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -44,10 +44,12 @@ def __init__( enable=True, path=None, dimensions=None, + block=None, ): super().__init__( sim_data, model_or_sim, structure, enable, path, dimensions ) + self._block = block if self.structure.layered: try: self._layer_shape = self.layer_shape() @@ -719,12 +721,18 @@ def _get_data(self, layer=None, apply_mult=False, **kwargs): layer = (layer,) storage = self._get_storage_obj() if storage is not None: + block_exists = self._block.header_exists( + self._current_key, self.path + ) try: - data = storage.get_data(layer, apply_mult) + data = storage.get_data( + layer, apply_mult, block_exists=block_exists + ) if ( "array" in kwargs and kwargs["array"] and isinstance(self, MFTransientArray) + and data != [] ): data = np.expand_dims(data, 0) return data @@ -1471,6 +1479,7 @@ def __init__( enable=True, path=None, dimensions=None, + block=None, ): super().__init__( sim_data=sim_data, @@ -1480,6 +1489,7 @@ def __init__( enable=enable, path=path, dimensions=dimensions, + block=block, ) self._transient_setup(self._data_storage) self.repeating = True @@ -1632,7 +1642,7 @@ def _get_array(self, num_sp, apply_mult, **kwargs): break if data is not None: if self.structure.type == DatumType.integer: - data = np.full_like(data, 0) + data = np.full_like(data, 1) else: data = np.full_like(data, 0.0) data = np.expand_dims(data, 0) @@ -1642,6 +1652,17 @@ def _get_array(self, num_sp, apply_mult, **kwargs): output = np.concatenate((output, data)) return output + def has_data(self, layer=None): + if layer is None: + for sto_key in self._data_storage.keys(): + self.get_data_prep(sto_key) + if super().has_data(): + return True + return False + else: + self.get_data_prep(layer) + return super().has_data() + def get_data(self, layer=None, apply_mult=True, **kwargs): """Returns the data associated with stress period key `layer`. If `layer` is None, returns all data for time `layer`. @@ -1677,6 +1698,10 @@ def get_data(self, layer=None, apply_mult=True, **kwargs): else: output = {sp: data} else: + if data is None and self._block.header_exists( + self._current_key, self.path + ): + data = [] if "array" in kwargs: output.append(data) else: diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index 400e157408..02cadf92a9 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -49,6 +49,7 @@ def __init__( path=None, dimensions=None, package=None, + block=None, ): super().__init__( sim_data, model_or_sim, structure, enable, path, dimensions @@ -72,6 +73,7 @@ def __init__( ex, ) self._package = package + self._block = block self._last_line_info = [] self._data_line = None self._temp_dict = {} @@ -316,7 +318,7 @@ def store_internal( } self._set_data(internal_data, check_data=check_data) - def has_data(self): + def has_data(self, key=None): """Returns whether this MFList has any data associated with it.""" try: if self._get_storage_obj() is None: @@ -343,7 +345,10 @@ def _get_data(self, apply_mult=False, **kwargs): try: if self._get_storage_obj() is None: return None - return self._get_storage_obj().get_data() + block_exists = self._block.header_exists( + self._current_key, self.path + ) + return self._get_storage_obj().get_data(block_exists=block_exists) except Exception as ex: type_, value_, traceback_ = sys.exc_info() raise MFDataException( @@ -1393,6 +1398,7 @@ def __init__( path=None, dimensions=None, package=None, + block=None, ): super().__init__( sim_data=sim_data, @@ -1403,6 +1409,7 @@ def __init__( path=path, dimensions=dimensions, package=package, + block=block, ) self._transient_setup(self._data_storage) self.repeating = True @@ -1618,6 +1625,19 @@ def store_internal( ) self._cache_model_grid = False + def has_data(self, key=None): + """Returns whether this MFList has any data associated with it in key + "key".""" + if key is None: + for sto_key in self._data_storage.keys(): + self.get_data_prep(sto_key) + if super().has_data(): + return True + return False + else: + self.get_data_prep(layer) + return super().has_data() + def get_data(self, key=None, apply_mult=False, **kwargs): """Returns the data for stress period `key`. @@ -1646,6 +1666,8 @@ def get_data(self, key=None, apply_mult=False, **kwargs): if sp in self._data_storage: self.get_data_prep(sp) data = super().get_data(apply_mult=apply_mult) + elif self._block.header_exists(sp): + data = None output.append(data) return output else: @@ -1966,6 +1988,7 @@ def __init__( path=None, dimensions=None, package=None, + block=None, ): super().__init__( sim_data=sim_data, @@ -1975,6 +1998,7 @@ def __init__( path=path, dimensions=dimensions, package=package, + block=block, ) def get_data(self, key=None, apply_mult=False, **kwargs): diff --git a/flopy/mf6/data/mfdatascalar.py b/flopy/mf6/data/mfdatascalar.py index 6bab86163b..11e1327369 100644 --- a/flopy/mf6/data/mfdatascalar.py +++ b/flopy/mf6/data/mfdatascalar.py @@ -82,7 +82,7 @@ def dtype(self): return np.int32 return None - def has_data(self): + def has_data(self, key=None): """Returns whether this object has data associated with it.""" try: return self._get_storage_obj().has_data() @@ -771,16 +771,13 @@ def add_one(self, key=0): def has_data(self, key=None): if key is None: - data_found = False for sto_key in self._data_storage.keys(): self.get_data_prep(sto_key) - data_found = data_found or super().has_data() - if data_found: - break + if super().has_data(): + return True else: self.get_data_prep(key) - data_found = super().has_data() - return data_found + return super().has_data() def get_data(self, key=0, **kwargs): """Returns the data for stress period `key`. @@ -852,7 +849,7 @@ def get_file_entry( ext_file_action=ext_file_action ) file_entry.append(text_entry) - if file_entry > 1: + if len(file_entry) > 1: return "\n\n".join(file_entry) elif file_entry == 1: return file_entry[0] diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py index c59cb3de5b..34f08be33b 100644 --- a/flopy/mf6/data/mfdatastorage.py +++ b/flopy/mf6/data/mfdatastorage.py @@ -633,10 +633,14 @@ def get_const_val(self, layer=None): def has_data(self, layer=None): ret_val = self._access_data(layer, False) - return ret_val is not None and ret_val != False + return ret_val is not None and ret_val is not False - def get_data(self, layer=None, apply_mult=True): - return self._access_data(layer, True, apply_mult=apply_mult) + def get_data(self, layer=None, apply_mult=True, block_exists=False): + data = self._access_data(layer, True, apply_mult=apply_mult) + if data is None and block_exists: + return [] + else: + return data def _access_data(self, layer, return_data=False, apply_mult=True): layer_check = self._resolve_layer(layer) @@ -2341,7 +2345,10 @@ def _build_full_data(self, apply_multiplier=False): ): full_data = data_out else: - full_data[layer] = data_out + if is_aux and full_data.shape == data_out.shape: + full_data = data_out + else: + full_data[layer] = data_out if is_aux: if full_data is not None: all_none = False diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index 8b6543c551..2b2bce5325 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -56,10 +56,17 @@ class MFBlockHeader: """ def __init__( - self, name, variable_strings, comment, simulation_data=None, path=None + self, + name, + variable_strings, + comment, + simulation_data=None, + path=None, + block=None, ): self.name = name self.variable_strings = variable_strings + self.block = block if not ( (simulation_data is None and path is None) or (simulation_data is not None and path is not None) @@ -103,6 +110,17 @@ def __init__( self.blk_trailing_comment_path = ("blk_trailing_comment",) self.blk_post_comment_path = ("blk_post_comment",) + def __lt__(self, other): + transient_key = self.get_transient_key() + if transient_key is None: + return True + else: + other_key = other.get_transient_key() + if other_key is None: + return False + else: + return transient_key < other_key + def build_header_variables( self, simulation_data, @@ -132,7 +150,7 @@ def build_header_variables( if len(fixed_data) > 0: fixed_data = [tuple(fixed_data)] # create data object - new_data = MFBlock.data_factory( + new_data = self.block.data_factory( simulation_data, None, block_header_structure[0], @@ -257,11 +275,14 @@ def write_footer(self, fd): fd.write(str(entry.rstrip())) fd.write("\n") - def get_transient_key(self): + def get_transient_key(self, data_path=None): """Get transient key associated with this block header.""" transient_key = None for index in range(0, len(self.data_items)): if self.data_items[index].structure.type != DatumType.keyword: + if data_path == self.data_items[index].path: + # avoid infinite recursion + return True transient_key = self.data_items[index].get_data() if isinstance(transient_key, np.recarray): item_struct = self.data_items[index].structure @@ -342,6 +363,7 @@ def __init__( MFComment("", path, simulation_data, 0), simulation_data, path, + self, ) ] self.structure = structure @@ -378,8 +400,8 @@ def _get_data_str(self, formal): return data_str # return an MFScalar, MFList, or MFArray - @staticmethod def data_factory( + self, sim_data, model_or_sim, structure, @@ -424,10 +446,17 @@ def data_factory( enable, path, dimensions, + self, ) elif data_type == mfstructure.DataType.array_transient: trans_array = mfdataarray.MFTransientArray( - sim_data, model_or_sim, structure, enable, path, dimensions + sim_data, + model_or_sim, + structure, + enable, + path, + dimensions, + self, ) if data is not None: trans_array.set_data(data, key=0) @@ -442,6 +471,7 @@ def data_factory( path, dimensions, package, + self, ) elif data_type == mfstructure.DataType.list_transient: trans_list = mfdatalist.MFTransientList( @@ -452,6 +482,7 @@ def data_factory( path, dimensions, package, + self, ) if data is not None: trans_list.set_data(data, key=0, autofill=True) @@ -465,6 +496,7 @@ def data_factory( path, dimensions, package, + self, ) if data is not None: mult_list.set_data(data, key=0, autofill=True) @@ -594,7 +626,7 @@ def add_dataset(self, dataset_struct, data, var_path): return self.datasets[var_path[-1]] def _build_repeating_header(self, header_data): - if self._header_exists(header_data[0]): + if self.header_exists(header_data[0]): return if ( len(self.block_headers[-1].data_items) == 1 @@ -607,6 +639,7 @@ def _build_repeating_header(self, header_data): MFComment("", self.path, self._simulation_data, 0), self._simulation_data, block_header_path, + self, ) self.block_headers.append(block_header) else: @@ -651,7 +684,7 @@ def _new_dataset( else: initial_val_path = initial_val try: - new_data = MFBlock.data_factory( + new_data = self.data_factory( self._simulation_data, self._model_or_sim, dataset_struct, @@ -765,6 +798,7 @@ def load(self, block_header, fd, strict=True): self.enabled = True if not self.loaded: self.block_headers = [] + block_header.block = self self.block_headers.append(block_header) # process any header variable @@ -1154,25 +1188,27 @@ def write(self, fd, ext_file_action=ExtFileAction.copy_relative_paths): for repeating_dataset in repeating_datasets: # resolve any missing block headers self._add_missing_block_headers(repeating_dataset) - for block_header in self.block_headers: + for block_header in sorted(self.block_headers): # write block self._write_block(fd, block_header, ext_file_action) - else: self._write_block(fd, self.block_headers[0], ext_file_action) def _add_missing_block_headers(self, repeating_dataset): for key in repeating_dataset.get_active_key_list(): - if not self._header_exists(key[0]): + data = repeating_dataset.get_data(key) + if not self.header_exists(key[0]) and data is not None: self._build_repeating_header([key[0]]) - def _header_exists(self, key): + def header_exists(self, key, data_path=None): if not isinstance(key, list): comp_key_list = [key] else: comp_key_list = key for block_header in self.block_headers: - transient_key = block_header.get_transient_key() + transient_key = block_header.get_transient_key(data_path) + if transient_key is True: + return for comp_key in comp_key_list: if transient_key is not None and transient_key == comp_key: return True @@ -1259,40 +1295,13 @@ def _find_repeating_datasets(self): return repeating_datasets def _write_block(self, fd, block_header, ext_file_action): - # write block header - block_header.write_header(fd) transient_key = None if len(block_header.data_items) > 0: transient_key = block_header.get_transient_key() - if self.external_file_name is not None: - # write block contents to external file - indent_string = self._simulation_data.indent_string - fd.write(f"{indent_string}open/close {self.external_file_name}\n") - fd_main = fd - fd_path = os.path.split(os.path.realpath(fd.name))[0] - try: - fd = open(os.path.join(fd_path, self.external_file_name), "w") - except: - type_, value_, traceback_ = sys.exc_info() - message = ( - f'Error reading external file "{self.external_file_name}"' - ) - raise MFDataException( - self._container_package.model_name, - self._container_package._get_pname(), - self.path, - "reading external file", - self.structure.name, - inspect.stack()[0][3], - type_, - value_, - traceback_, - message, - self._simulation_data.debug, - ) - - # write data sets + # gather data sets to write + data_set_output = [] + data_found = False for key, dataset in self.datasets.items(): try: if transient_key is None: @@ -1303,9 +1312,10 @@ def _write_block(self, fd, block_header, ext_file_action): print( f" writing data {dataset.structure.name}..." ) - fd.write( + data_set_output.append( dataset.get_file_entry(ext_file_action=ext_file_action) ) + data_found = True else: if ( self._simulation_data.verbosity_level.value @@ -1316,17 +1326,19 @@ def _write_block(self, fd, block_header, ext_file_action): ".".format(dataset.structure.name, transient_key) ) if dataset.repeating: - fd.write( - dataset.get_file_entry( - transient_key, ext_file_action=ext_file_action - ) + output = dataset.get_file_entry( + transient_key, ext_file_action=ext_file_action ) + if output is not None: + data_set_output.append(output) + data_found = True else: - fd.write( + data_set_output.append( dataset.get_file_entry( ext_file_action=ext_file_action ) ) + data_found = True except MFDataException as mfde: raise MFDataException( mfdata_except=mfde, @@ -1338,6 +1350,42 @@ def _write_block(self, fd, block_header, ext_file_action): f'"{self.structure.name}" to file "{fd.name}"' ), ) + if not data_found: + return + # write block header + block_header.write_header(fd) + + if self.external_file_name is not None: + # write block contents to external file + indent_string = self._simulation_data.indent_string + fd.write(f"{indent_string}open/close {self.external_file_name}\n") + fd_main = fd + fd_path = os.path.split(os.path.realpath(fd.name))[0] + try: + fd = open(os.path.join(fd_path, self.external_file_name), "w") + except: + type_, value_, traceback_ = sys.exc_info() + message = ( + f'Error reading external file "{self.external_file_name}"' + ) + raise MFDataException( + self._container_package.model_name, + self._container_package._get_pname(), + self.path, + "reading external file", + self.structure.name, + inspect.stack()[0][3], + type_, + value_, + traceback_, + message, + self._simulation_data.debug, + ) + + # write data sets + for output in data_set_output: + fd.write(output) + # write trailing comments pth = block_header.blk_trailing_comment_path if pth in self._simulation_data.mfdata: From 50e51bf3c14bbfea78a480399ffe36fdfbfd2451 Mon Sep 17 00:00:00 2001 From: Scott Paulinski Date: Mon, 29 Nov 2021 12:57:05 -0800 Subject: [PATCH 4/6] fix(array and getdata): added tests for array and getdata for different types of empty stress periods --- autotest/t505_test.py | 179 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 167 insertions(+), 12 deletions(-) diff --git a/autotest/t505_test.py b/autotest/t505_test.py index 81193ed875..3679025ab2 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -353,12 +353,19 @@ def test_multi_model(): def test_array(): + # get_data + # empty data in period block vs data repeating + # array + # aux values, test that they work the same as other arrays (is a value + # of zero always used even if aux is defined in a previous stress + # period?) + import flopy mf6 = flopy.mf6 - sim_name = "testsim" - model_name = "testmodel" - out_dir = "tmp" + sim_name = "test_array" + model_name = "test_array" + out_dir = os.path.join("temp", "t505_test_array") tdis_name = "{}.tdis".format(sim_name) sim = mf6.MFSimulation( @@ -366,7 +373,22 @@ def test_array(): ) tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0), (6.0, 3, 1.0), (6.0, 3, 1.0)] tdis = mf6.ModflowTdis(sim, time_units="DAYS", nper=4, perioddata=tdis_rc) - + ims_package = ModflowIms( + sim, + pname="my_ims_file", + filename=f"{sim_name}.ims", + print_option="ALL", + complexity="SIMPLE", + outer_dvclose=0.00001, + outer_maximum=50, + under_relaxation="NONE", + inner_maximum=30, + inner_dvclose=0.00001, + linear_acceleration="CG", + preconditioner_levels=7, + preconditioner_drop_tolerance=0.01, + number_orthogonalizations=2, + ) model = mf6.ModflowGwf( sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) ) @@ -374,27 +396,66 @@ def test_array(): dis = mf6.ModflowGwfdis( model, length_units="FEET", - nlay=1, + nlay=4, nrow=2, ncol=2, delr=500.0, delc=500.0, top=100.0, - botm=50.0, + botm=[50.0, 0.0, -50.0, -100.0], filename="{}.dis".format(model_name), ) - irch = {1: [[0, 1], [2, 1]], 2: [[0, 1], [2, 3]]} - rcha = mf6.ModflowGwfrcha(model, irch=irch, recharge={0: 1.0, 2: 2.0}) + ic_package = mf6.ModflowGwfic( + model, strt=90.0, filename=f"{model_name}.ic" + ) + npf_package = mf6.ModflowGwfnpf( + model, + pname="npf_1", + save_flows=True, + alternative_cell_averaging="logarithmic", + icelltype=1, + k=5.0, + ) + + aux = {1: [[50.0], [1.3]], 3: [[200.0], [1.5]]} + irch = {1: [[0, 2], [2, 1]], 2: [[0, 1], [2, 3]]} + rcha = mf6.ModflowGwfrcha( + model, + print_input=True, + print_flows=True, + auxiliary=[("var1", "var2")], + irch=irch, + recharge={1: 1.0, 2: 2.0}, + aux=aux, + ) val_irch = rcha.irch.array.sum(axis=(1, 2, 3)) - assert val_irch[0] == 0 - assert val_irch[1] == 4 + assert val_irch[0] == 4 + assert val_irch[1] == 5 assert val_irch[2] == 6 assert val_irch[3] == 6 + val_irch_2 = rcha.irch.get_data() + assert val_irch_2[0] is None + assert val_irch_2[1][1, 1] == 1 + assert val_irch_2[2][1, 1] == 3 + assert val_irch_2[3] == [] val_rch = rcha.recharge.array.sum(axis=(1, 2, 3)) - assert val_rch[0] == 4.0 + assert val_rch[0] == 0.0 assert val_rch[1] == 4.0 assert val_rch[2] == 8.0 assert val_rch[3] == 8.0 + val_rch_2 = rcha.recharge.get_data() + assert val_rch_2[0] is None + assert val_rch_2[1][0, 0] == 1.0 + assert val_rch_2[2][0, 0] == 2.0 + assert val_rch_2[3] == [] + aux_data_0 = rcha.aux.get_data(0) + assert aux_data_0 is None + aux_data_1 = rcha.aux.get_data(1) + assert aux_data_1[0][0][0] == 50.0 + aux_data_2 = rcha.aux.get_data(2) + assert aux_data_2 == [] + aux_data_3 = rcha.aux.get_data(3) + assert aux_data_3[0][0][0] == 200.0 welspdict = {1: [[(0, 0, 0), -25.0, 0.0]], 2: [[(0, 0, 0), 25.0, 0.0]]} wel = flopy.mf6.ModflowGwfwel( @@ -412,6 +473,99 @@ def test_array(): assert wel_array[2][0][1] == 25.0 assert wel_array[3][0][1] == 25.0 + drnspdict = { + 0: [[(0, 0, 0), 60.0, 10.0]], + 2: [], + 3: [[(0, 0, 0), 55.0, 5.0]], + } + drn = flopy.mf6.ModflowGwfdrn( + model, + print_input=True, + print_flows=True, + stress_period_data=drnspdict, + save_flows=False, + pname="DRN-1", + ) + drn_array = drn.stress_period_data.array + assert drn_array[0][0][1] == 60.0 + assert drn_array[1][0][1] == 60.0 + assert drn_array[2] is None + assert drn_array[3][0][1] == 55.0 + drn_gd_0 = drn.stress_period_data.get_data(0) + assert drn_gd_0[0][1] == 60.0 + drn_gd_1 = drn.stress_period_data.get_data(1) + assert drn_gd_1 is None + drn_gd_2 = drn.stress_period_data.get_data(2) + assert drn_gd_2 == [] + drn_gd_3 = drn.stress_period_data.get_data(3) + assert drn_gd_3[0][1] == 55.0 + + # test writing and loading model + sim.write_simulation() + if run: + sim.run_simulation() + + test_sim = MFSimulation.load( + sim_name, + "mf6", + exe_name, + out_dir, + write_headers=False, + ) + model = test_sim.get_model() + rcha = model.get_package("rcha") + wel = model.get_package("wel") + drn = model.get_package("drn") + # do same tests as above + val_irch = rcha.irch.array.sum(axis=(1, 2, 3)) + assert val_irch[0] == 4 + assert val_irch[1] == 5 + assert val_irch[2] == 6 + assert val_irch[3] == 6 + val_irch_2 = rcha.irch.get_data() + assert val_irch_2[0] is None + assert val_irch_2[1][1, 1] == 1 + assert val_irch_2[2][1, 1] == 3 + assert val_irch_2[3] == [] + val_rch = rcha.recharge.array.sum(axis=(1, 2, 3)) + assert val_rch[0] == 0.0 + assert val_rch[1] == 4.0 + assert val_rch[2] == 8.0 + assert val_rch[3] == 8.0 + val_rch_2 = rcha.recharge.get_data() + assert val_rch_2[0] is None + assert val_rch_2[1][0, 0] == 1.0 + assert val_rch_2[2][0, 0] == 2.0 + assert val_rch_2[3] == [] + aux_data_0 = rcha.aux.get_data(0) + assert aux_data_0 is None + aux_data_1 = rcha.aux.get_data(1) + assert aux_data_1[0][0][0] == 50.0 + aux_data_2 = rcha.aux.get_data(2) + assert aux_data_2 == [] + aux_data_3 = rcha.aux.get_data(3) + assert aux_data_3[0][0][0] == 200.0 + + wel_array = wel.stress_period_data.array + assert wel_array[0] is None + assert wel_array[1][0][1] == -25.0 + assert wel_array[2][0][1] == 25.0 + assert wel_array[3][0][1] == 25.0 + + drn_array = drn.stress_period_data.array + assert drn_array[0][0][1] == 60.0 + assert drn_array[1][0][1] == 60.0 + assert drn_array[2] is None + assert drn_array[3][0][1] == 55.0 + drn_gd_0 = drn.stress_period_data.get_data(0) + assert drn_gd_0[0][1] == 60.0 + drn_gd_1 = drn.stress_period_data.get_data(1) + assert drn_gd_1 is None + drn_gd_2 = drn.stress_period_data.get_data(2) + assert drn_gd_2 == [] + drn_gd_3 = drn.stress_period_data.get_data(3) + assert drn_gd_3[0][1] == 55.0 + def test_np001(): # init paths @@ -819,6 +973,7 @@ def test_np001(): sim.rename_all_packages("file_rename") sim.set_sim_path(rename_folder) sim.write_simulation() + if run: sim.run_simulation() sim.delete_output_files() @@ -1222,7 +1377,7 @@ def test_np002(): md2 = sim2.get_model() ghb2 = md2.get_package("ghb") spd2 = ghb2.stress_period_data.get_data(1) - assert spd2 is None + assert spd2 == [] # test paths sim_path_test = os.path.join(run_folder, "sim_path") From 1bde6008bd4ad4c8b58e158df7070e442d78eb3d Mon Sep 17 00:00:00 2001 From: Scott Paulinski Date: Tue, 30 Nov 2021 08:01:22 -0800 Subject: [PATCH 5/6] fix(array output): minor fix for how arrays output data --- flopy/mf6/data/mfdataarray.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index 0ff3fbc9f8..cd4a24fa7e 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -721,18 +721,13 @@ def _get_data(self, layer=None, apply_mult=False, **kwargs): layer = (layer,) storage = self._get_storage_obj() if storage is not None: - block_exists = self._block.header_exists( - self._current_key, self.path - ) try: - data = storage.get_data( - layer, apply_mult, block_exists=block_exists - ) + data = storage.get_data(layer, apply_mult) if ( "array" in kwargs and kwargs["array"] and isinstance(self, MFTransientArray) - and data != [] + and data is not [] ): data = np.expand_dims(data, 0) return data @@ -1698,10 +1693,6 @@ def get_data(self, layer=None, apply_mult=True, **kwargs): else: output = {sp: data} else: - if data is None and self._block.header_exists( - self._current_key, self.path - ): - data = [] if "array" in kwargs: output.append(data) else: From 6d2a6d2e8df5f592fa7406552b26ee075fc1db8d Mon Sep 17 00:00:00 2001 From: Scott Paulinski Date: Tue, 30 Nov 2021 08:02:47 -0800 Subject: [PATCH 6/6] fix(test): fixed convergence issues with new test case --- autotest/t505_test.py | 67 ++++++++++++++++++++++++++----------------- 1 file changed, 40 insertions(+), 27 deletions(-) diff --git a/autotest/t505_test.py b/autotest/t505_test.py index 3679025ab2..71953c946d 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -399,8 +399,8 @@ def test_array(): nlay=4, nrow=2, ncol=2, - delr=500.0, - delc=500.0, + delr=5000.0, + delc=5000.0, top=100.0, botm=[50.0, 0.0, -50.0, -100.0], filename="{}.dis".format(model_name), @@ -414,7 +414,7 @@ def test_array(): save_flows=True, alternative_cell_averaging="logarithmic", icelltype=1, - k=5.0, + k=50.0, ) aux = {1: [[50.0], [1.3]], 3: [[200.0], [1.5]]} @@ -425,7 +425,7 @@ def test_array(): print_flows=True, auxiliary=[("var1", "var2")], irch=irch, - recharge={1: 1.0, 2: 2.0}, + recharge={1: 0.0001, 2: 0.0002}, aux=aux, ) val_irch = rcha.irch.array.sum(axis=(1, 2, 3)) @@ -437,27 +437,29 @@ def test_array(): assert val_irch_2[0] is None assert val_irch_2[1][1, 1] == 1 assert val_irch_2[2][1, 1] == 3 - assert val_irch_2[3] == [] + assert val_irch_2[3] is None + val_irch_2_3 = rcha.irch.get_data(3) + assert val_irch_2_3 is None val_rch = rcha.recharge.array.sum(axis=(1, 2, 3)) assert val_rch[0] == 0.0 - assert val_rch[1] == 4.0 - assert val_rch[2] == 8.0 - assert val_rch[3] == 8.0 + assert val_rch[1] == 0.0004 + assert val_rch[2] == 0.0008 + assert val_rch[3] == 0.0008 val_rch_2 = rcha.recharge.get_data() assert val_rch_2[0] is None - assert val_rch_2[1][0, 0] == 1.0 - assert val_rch_2[2][0, 0] == 2.0 - assert val_rch_2[3] == [] + assert val_rch_2[1][0, 0] == 0.0001 + assert val_rch_2[2][0, 0] == 0.0002 + assert val_rch_2[3] is None aux_data_0 = rcha.aux.get_data(0) assert aux_data_0 is None aux_data_1 = rcha.aux.get_data(1) assert aux_data_1[0][0][0] == 50.0 aux_data_2 = rcha.aux.get_data(2) - assert aux_data_2 == [] + assert aux_data_2 is None aux_data_3 = rcha.aux.get_data(3) assert aux_data_3[0][0][0] == 200.0 - welspdict = {1: [[(0, 0, 0), -25.0, 0.0]], 2: [[(0, 0, 0), 25.0, 0.0]]} + welspdict = {1: [[(0, 0, 0), 0.25, 0.0]], 2: [[(0, 0, 0), 0.1, 0.0]]} wel = flopy.mf6.ModflowGwfwel( model, print_input=True, @@ -469,9 +471,9 @@ def test_array(): ) wel_array = wel.stress_period_data.array assert wel_array[0] is None - assert wel_array[1][0][1] == -25.0 - assert wel_array[2][0][1] == 25.0 - assert wel_array[3][0][1] == 25.0 + assert wel_array[1][0][1] == 0.25 + assert wel_array[2][0][1] == 0.1 + assert wel_array[3][0][1] == 0.1 drnspdict = { 0: [[(0, 0, 0), 60.0, 10.0]], @@ -500,6 +502,17 @@ def test_array(): drn_gd_3 = drn.stress_period_data.get_data(3) assert drn_gd_3[0][1] == 55.0 + ghbspdict = { + 0: [[(0, 1, 1), 60.0, 10.0]], + } + ghb = flopy.mf6.ModflowGwfghb( + model, + print_input=True, + print_flows=True, + stress_period_data=ghbspdict, + save_flows=False, + pname="GHB-1", + ) # test writing and loading model sim.write_simulation() if run: @@ -526,31 +539,31 @@ def test_array(): assert val_irch_2[0] is None assert val_irch_2[1][1, 1] == 1 assert val_irch_2[2][1, 1] == 3 - assert val_irch_2[3] == [] + assert val_irch_2[3] is None val_rch = rcha.recharge.array.sum(axis=(1, 2, 3)) assert val_rch[0] == 0.0 - assert val_rch[1] == 4.0 - assert val_rch[2] == 8.0 - assert val_rch[3] == 8.0 + assert val_rch[1] == 0.0004 + assert val_rch[2] == 0.0008 + assert val_rch[3] == 0.0008 val_rch_2 = rcha.recharge.get_data() assert val_rch_2[0] is None - assert val_rch_2[1][0, 0] == 1.0 - assert val_rch_2[2][0, 0] == 2.0 - assert val_rch_2[3] == [] + assert val_rch_2[1][0, 0] == 0.0001 + assert val_rch_2[2][0, 0] == 0.0002 + assert val_rch_2[3] is None aux_data_0 = rcha.aux.get_data(0) assert aux_data_0 is None aux_data_1 = rcha.aux.get_data(1) assert aux_data_1[0][0][0] == 50.0 aux_data_2 = rcha.aux.get_data(2) - assert aux_data_2 == [] + assert aux_data_2 is None aux_data_3 = rcha.aux.get_data(3) assert aux_data_3[0][0][0] == 200.0 wel_array = wel.stress_period_data.array assert wel_array[0] is None - assert wel_array[1][0][1] == -25.0 - assert wel_array[2][0][1] == 25.0 - assert wel_array[3][0][1] == 25.0 + assert wel_array[1][0][1] == 0.25 + assert wel_array[2][0][1] == 0.1 + assert wel_array[3][0][1] == 0.1 drn_array = drn.stress_period_data.array assert drn_array[0][0][1] == 60.0