From cc99fcd4b9315466f01ae7d6108160b5519d3cf6 Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Fri, 21 Aug 2020 13:25:45 -0700 Subject: [PATCH 1/6] feat(cellid not as tuple): list data input may now have cellids as separate integers rather than as integers in a tuple. flopy still internally stores the integers as a tuple. --- flopy/mf6/data/mfdatalist.py | 2 +- flopy/mf6/data/mfdatastorage.py | 185 ++++++++++++++++++++++---------- flopy/mf6/data/mfdatautil.py | 60 +++++++++-- 3 files changed, 181 insertions(+), 66 deletions(-) diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index b94bd869b2..9a0a31f57c 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -419,7 +419,7 @@ def _check_valid_cellids(self): data = storage_obj.get_data() # check data for invalid cellids for index, is_cellid in enumerate( - storage_obj.recarray_cellid_list + storage_obj.resolve_cellidlist(data) ): if is_cellid: for record in data: diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py index afaa314eb1..cf67ade0bb 100644 --- a/flopy/mf6/data/mfdatastorage.py +++ b/flopy/mf6/data/mfdatastorage.py @@ -730,9 +730,8 @@ def _access_data(self, layer, return_data=False, apply_mult=True): ): data_line = data_line + (None,) data_list.append(data_line) - return np.rec.array( - data_list, self._recarray_type_list - ) + type_list = self.resolve_typelist(data_list) + return np.rec.array(data_list, type_list) else: return ( self.layer_storage[ @@ -801,6 +800,7 @@ def append_data(self, data): self._simulation_data.debug, ) internal_data = self.layer_storage.first_item().internal_data + if internal_data is None: if len(data[0]) != len(self._recarray_type_list): # rebuild type list using existing data as a guide @@ -1187,6 +1187,10 @@ def store_internal( def _build_recarray(self, data, key, autofill): self.build_type_list(data=data, key=key) + if not self.tuple_cellids(data): + # fix data so cellid is a single tuple + data = self.make_tuple_cellids(data) + if autofill and data is not None: # resolve any fields with data types that do not # agree with the expected type list @@ -1236,6 +1240,54 @@ def _build_recarray(self, data, key, autofill): self._verify_list(new_data) return new_data + def make_tuple_cellids(self, data): + # convert cellids from individual layer, row, column fields into + # tuples (layer, row, column) + data_dim = self.data_dimensions + model_grid = data_dim.get_model_grid() + cellid_size = model_grid.get_num_spatial_coordinates() + + new_data = [] + current_cellid = () + for index, line in enumerate(data): + new_line = [] + for item, is_cellid in zip(line, self.recarray_cellid_list_ex): + if is_cellid: + current_cellid += (item,) + if len(current_cellid) == cellid_size: + new_line.append(current_cellid) + current_cellid = () + else: + new_line.append(item) + new_data.append(tuple(new_line)) + return new_data + + def tuple_cellids(self, data): + for data_entry, cellid in zip(data[0], self.recarray_cellid_list): + if cellid: + if isinstance(data_entry, int): + # cellid is stored in separate columns in the recarray + # (eg: one column for layer one column for row and + # one columne for column) + return False + else: + # cellid is stored in a single column in the recarray + # as a tuple + return True + return True + + def resolve_typelist(self, data): + if self.tuple_cellids(data): + return self._recarray_type_list + else: + return self._recarray_type_list_ex + + def resolve_cellidlist(self, data): + if self.tuple_cellids(data): + return self.recarray_cellid_list + else: + return self.recarray_cellid_list_ex + def _resolve_multitype_fields(self, data): # find any data fields where the data is not a consistent type itype_len = len(self._recarray_type_list) @@ -1281,6 +1333,17 @@ def store_external( fp = self._simulation_data.mfpath.resolve_path( file_path, model_name ) + # store data internally first so that a file entry + # can be generated + self.store_internal( + data, + layer_new, + False, + [multiplier], + None, + False, + print_format, + ) if binary: file_access = MFFileAccessList( self.data_dimensions.structure, @@ -1290,7 +1353,7 @@ def store_external( self._stress_period, ) file_access.write_binary_file( - data, + self.layer_storage.first_item().internal_data, fp, self._model_or_sim.modeldiscrit, precision="double", @@ -1318,18 +1381,6 @@ def store_external( message, self._simulation_data.debug, ) - # store data internally first so that a file entry - # can be generated - self.store_internal( - data, - layer_new, - False, - [multiplier], - None, - False, - print_format, - ) - ext_file_entry = self._get_file_entry() fd.write(ext_file_entry) fd.close() @@ -2313,10 +2364,14 @@ def build_type_list( resolve_data_shape=True, key=None, nseg=None, + cellid_expanded=False, ): if data_set is None: + self.jagged_record = False self._recarray_type_list = [] + self._recarray_type_list_ex = [] self.recarray_cellid_list = [] + self.recarray_cellid_list_ex = [] data_set = self.data_dimensions.structure initial_keyword = True package_dim = self.data_dimensions.package_dim @@ -2343,50 +2398,42 @@ def build_type_list( if aux_var_names is not None: for aux_var_name in aux_var_names[0]: if aux_var_name.lower() != "auxiliary": - self._recarray_type_list.append( - (aux_var_name, data_type) + self._append_type_lists( + aux_var_name, data_type, False ) - self.recarray_cellid_list.append(False) elif data_item.type == DatumType.record: # record within a record, recurse self.build_type_list(data_item, True, data) elif data_item.type == DatumType.keystring: - self._recarray_type_list.append( - (data_item.name, data_type) + self.jagged_record = True + self._append_type_lists( + data_item.name, data_type, data_item.is_cellid ) - self.recarray_cellid_list.append(data_item.is_cellid) # add potential data after keystring to type list ks_data_item = deepcopy(data_item) ks_data_item.type = DatumType.string ks_data_item.name = "{}_data".format(ks_data_item.name) ks_rec_type = ks_data_item.get_rec_type() - self._recarray_type_list.append( - (ks_data_item.name, ks_rec_type) + self._append_type_lists( + ks_data_item.name, ks_rec_type, ks_data_item.is_cellid ) - self.recarray_cellid_list.append(ks_data_item.is_cellid) - if index == len(data_set.data_item_structures) - 1: + if ( + index == len(data_set.data_item_structures) - 1 + and data is not None + ): idx = 1 data_line_max_size = self._get_max_data_line_size(data) - while ( - data is not None - and len(self._recarray_type_list) - < data_line_max_size - ): + type_list = self.resolve_typelist(data) + while len(type_list) < data_line_max_size: # keystrings at the end of a line can contain items # of variable length. assume everything at the # end of the data line is related to the last # keystring - self._recarray_type_list.append( - ( - "{}_{}".format(ks_data_item.name, idx), - ks_rec_type, - ) + name = "{}_{}".format(ks_data_item.name, idx) + self._append_type_lists( + name, ks_rec_type, ks_data_item.is_cellid ) - self.recarray_cellid_list.append( - ks_data_item.is_cellid - ) - idx += 1 elif ( @@ -2406,11 +2453,9 @@ def build_type_list( data_item.type != DatumType.string and data_item.type != DatumType.keyword ): - self._recarray_type_list.append( - ("{}_label".format(data_item.name), object) - ) - self.recarray_cellid_list.append( - data_item.is_cellid + name = "{}_label".format(data_item.name) + self._append_type_lists( + name, object, data_item.is_cellid ) if ( nseg is not None @@ -2474,6 +2519,7 @@ def build_type_list( # shape is indeterminate 1-d array and no data # provided to resolve resolved_shape[0] = 1 + self.jagged_record = True if data_item.is_cellid: if ( data_item.shape is not None @@ -2488,22 +2534,16 @@ def build_type_list( data_item.remove_cellid(resolved_shape, size) for index in range(0, resolved_shape[0]): if resolved_shape[0] > 1: - # type list fields must have unique names - self._recarray_type_list.append( - ( - "{}_{}".format(data_item.name, index), - data_type, - ) - ) + name = "{}_{}".format(data_item.name, index) else: - self._recarray_type_list.append( - (data_item.name, data_type) - ) - self.recarray_cellid_list.append( - data_item.is_cellid + name = data_item.name + self._append_type_lists( + name, data_type, data_item.is_cellid ) - - return self._recarray_type_list + if cellid_expanded: + return self._recarray_type_list_ex + else: + return self._recarray_type_list def get_default_mult(self): if self._data_type == DatumType.integer: @@ -2511,6 +2551,33 @@ def get_default_mult(self): else: return 1.0 + def _append_type_lists(self, name, data_type, iscellid): + # add entry(s) to type lists + self._recarray_type_list.append((name, data_type)) + self.recarray_cellid_list.append(iscellid) + if iscellid and self._model_or_sim.model_type is not None: + # write each part of the cellid out as a separate entry + # to _recarray_list_list_ex + model_grid = self.data_dimensions.get_model_grid() + cellid_size = model_grid.get_num_spatial_coordinates() + # determine header for different grid types + if cellid_size == 1: + self._recarray_type_list_ex.append((name, int)) + elif cellid_size == 2: + self._recarray_type_list_ex.append(("layer", int)) + self._recarray_type_list_ex.append(("cell2d_num", int)) + else: + self._recarray_type_list_ex.append(("layer", int)) + self._recarray_type_list_ex.append(("row", int)) + self._recarray_type_list_ex.append(("column", int)) + + for index in range(0, cellid_size): + self.recarray_cellid_list_ex.append(iscellid) + else: + self._recarray_type_list_ex.append((name, data_type)) + self.recarray_cellid_list_ex.append(iscellid) + return iscellid + @staticmethod def _calc_data_size(data, count_to=None, current_length=None): if current_length is None: diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 530533c5d3..6aab1db991 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -678,20 +678,19 @@ def _build_template_data(self, type_list): template_data.append(None) return tuple(template_data) - def empty( + def dtype( self, model, - maxbound=None, aux_vars=None, boundnames=False, nseg=None, timeseries=False, - stress_periods=None, + cellid_expanded=False, ): - from ..data import mfdatastorage, mfstructure + from ..data import mfdatastorage + # get data storage data_struct, data_dimensions = self._get_data_dimensions(model) - data_type = data_struct.get_datatype() # build a temporary data storage object data_storage = mfdatastorage.DataStorage( model.simulation_data, @@ -703,7 +702,33 @@ def empty( ) # build type list - type_list = data_storage.build_type_list(nseg=nseg) + type_list = data_storage.build_type_list( + nseg=nseg, cellid_expanded=cellid_expanded + ) + if data_storage.jagged_record: + comment = ( + "Data dimensions can not be determined for " + "{}. Data structure may be jagged or may contain " + "a keystring. Data type information is therefore " + "dependant on the data and can not be retreived " + "prior to the data being loaded" + ".".format(data_storage.data_dimensions.structure.name) + ) + type_, value_, traceback_ = sys.exc_info() + + raise MFDataException( + data_struct.get_model(), + data_struct.get_package(), + data_struct.path, + "generating array template", + data_struct.name, + inspect.stack()[0][3], + type_, + value_, + traceback_, + comment, + model.simulation_data.debug, + ) if aux_vars is not None: if len(aux_vars) > 0 and ( isinstance(aux_vars[0], list) or isinstance(aux_vars[0], tuple) @@ -718,6 +743,29 @@ def empty( # fix type list to make all types objects for index, d_type in enumerate(type_list): type_list[index] = (d_type[0], object) + return type_list + + def empty( + self, + model, + maxbound=None, + aux_vars=None, + boundnames=False, + nseg=None, + timeseries=False, + stress_periods=None, + cellid_expanded=False, + ): + from ..data import mfstructure + + # get type list + type_list = self.dtype( + model, aux_vars, boundnames, nseg, timeseries, cellid_expanded, + ) + + # get data storage + data_struct, data_dimensions = self._get_data_dimensions(model) + data_type = data_struct.get_datatype() # build recarray template_data = self._build_template_data(type_list) From 4bac7157e5c50fe0615e110c92c5fe9fad023478 Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Fri, 21 Aug 2020 13:27:35 -0700 Subject: [PATCH 2/6] feat(cellid not as tuple): added tests for this feature --- autotest/t505_test.py | 3122 +++++++++++++++++++++++++++-------------- 1 file changed, 2062 insertions(+), 1060 deletions(-) diff --git a/autotest/t505_test.py b/autotest/t505_test.py index bb04d0f174..a90dc882a5 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -39,16 +39,16 @@ try: import pymake except: - print('could not import pymake') + print("could not import pymake") -exe_name = 'mf6' +exe_name = "mf6" v = flopy.which(exe_name) run = True if v is None: run = False -cpth = os.path.join('temp', 't505') +cpth = os.path.join("temp", "t505") # make the directory if it does not exist if not os.path.isdir(cpth): os.makedirs(cpth) @@ -56,169 +56,275 @@ def np001(): # init paths - test_ex_name = 'np001' - model_name = 'np001_mod' + test_ex_name = "np001" + model_name = "np001_mod" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'np001_mod.hds') - expected_cbc_file = os.path.join(expected_output_folder, 'np001_mod.cbc') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "np001_mod.hds") + expected_cbc_file = os.path.join(expected_output_folder, "np001_mod.cbc") # model tests - test_sim = MFSimulation(sim_name=test_ex_name, version='mf6', - exe_name=exe_name, sim_ws=run_folder, - continue_=True, memory_print_option='summary') + test_sim = MFSimulation( + sim_name=test_ex_name, + version="mf6", + exe_name=exe_name, + sim_ws=run_folder, + continue_=True, + memory_print_option="summary", + ) name = test_sim.name_file assert name.continue_.get_data() assert name.nocheck.get_data() is None - assert name.memory_print_option.get_data() == 'summary' + assert name.memory_print_option.get_data() == "summary" kwargs = {} - kwargs['bad_kwarg'] = 20 + kwargs["bad_kwarg"] = 20 try: ex = False - bad_model = ModflowGwf(test_sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name), - **kwargs) + bad_model = ModflowGwf( + test_sim, + modelname=model_name, + model_nam_file="{}.nam".format(model_name), + **kwargs + ) except FlopyException: ex = True - assert (ex == True) + assert ex == True kwargs = {} - kwargs['xul'] = 20.5 - good_model = ModflowGwf(test_sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name), - model_rel_path='model_folder', - **kwargs) + kwargs["xul"] = 20.5 + good_model = ModflowGwf( + test_sim, + modelname=model_name, + model_nam_file="{}.nam".format(model_name), + model_rel_path="model_folder", + **kwargs + ) # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth, write_headers=False) + sim = MFSimulation( + sim_name=test_ex_name, + version="mf6", + exe_name=exe_name, + sim_ws=pth, + write_headers=False, + ) tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=1, - perioddata=[(2.0, 1, 1.0)]) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=1, perioddata=[(2.0, 1, 1.0)] + ) # specifying the tdis package twice should remove the old tdis package - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=2, - perioddata=tdis_rc) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=2, perioddata=tdis_rc + ) # first ims file to be replaced - ims_package = ModflowIms(sim, pname='my_ims_file', filename='old_name.ims', - print_option='ALL', complexity='SIMPLE', - outer_hclose=0.00001, - outer_maximum=10, under_relaxation='NONE', - inner_maximum=10, - inner_hclose=0.001, linear_acceleration='CG', - preconditioner_levels=2, - preconditioner_drop_tolerance=0.00001, - number_orthogonalizations=5) + ims_package = ModflowIms( + sim, + pname="my_ims_file", + filename="old_name.ims", + print_option="ALL", + complexity="SIMPLE", + outer_hclose=0.00001, + outer_maximum=10, + under_relaxation="NONE", + inner_maximum=10, + inner_hclose=0.001, + linear_acceleration="CG", + preconditioner_levels=2, + preconditioner_drop_tolerance=0.00001, + number_orthogonalizations=5, + ) # replace with real ims file - ims_package = ModflowIms(sim, pname='my_ims_file', - filename='{}.ims'.format(test_ex_name), - print_option='ALL', complexity='SIMPLE', - outer_hclose=0.00001, - outer_maximum=50, under_relaxation='NONE', - inner_maximum=30, - inner_hclose=0.00001, linear_acceleration='CG', - preconditioner_levels=7, - preconditioner_drop_tolerance=0.01, - number_orthogonalizations=2) - - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) + ims_package = ModflowIms( + sim, + pname="my_ims_file", + filename="{}.ims".format(test_ex_name), + print_option="ALL", + complexity="SIMPLE", + outer_hclose=0.00001, + outer_maximum=50, + under_relaxation="NONE", + inner_maximum=30, + inner_hclose=0.00001, + linear_acceleration="CG", + preconditioner_levels=7, + preconditioner_drop_tolerance=0.01, + number_orthogonalizations=2, + ) + + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) # test case insensitive lookup - assert(sim.get_model(model_name.upper()) is not None) + assert sim.get_model(model_name.upper()) is not None # test getting model using attribute model = sim.np001_mod - assert(model is not None and model.name == 'np001_mod') + assert model is not None and model.name == "np001_mod" tdis = sim.tdis - assert(tdis is not None and tdis.package_type == 'tdis') - - dis_package = flopy.mf6.ModflowGwfdis(model, length_units='FEET', nlay=1, - nrow=1, ncol=1, delr=100.0, - delc=100.0, - top=60.0, botm=50.0, - filename='{}.dis'.format(model_name), - pname='mydispkg') + assert tdis is not None and tdis.package_type == "tdis" + + dis_package = flopy.mf6.ModflowGwfdis( + model, + length_units="FEET", + nlay=1, + nrow=1, + ncol=1, + delr=100.0, + delc=100.0, + top=60.0, + botm=50.0, + filename="{}.dis".format(model_name), + pname="mydispkg", + ) # specifying dis package twice with the same name should automatically # remove the old dis package - top = {'filename': 'top.bin', 'data': 100.0, 'binary': True} - botm = {'filename': 'botm.bin', 'data': 50.0, 'binary': True} - dis_package = flopy.mf6.ModflowGwfdis(model, length_units='FEET', nlay=1, - nrow=1, ncol=10, delr=500.0, - delc=500.0, - top=top, botm=botm, - filename='{}.dis'.format(model_name), - pname='mydispkg') + top = {"filename": "top.bin", "data": 100.0, "binary": True} + botm = {"filename": "botm.bin", "data": 50.0, "binary": True} + dis_package = flopy.mf6.ModflowGwfdis( + model, + length_units="FEET", + nlay=1, + nrow=1, + ncol=10, + delr=500.0, + delc=500.0, + top=top, + botm=botm, + filename="{}.dis".format(model_name), + pname="mydispkg", + ) top_data = dis_package.top.get_data() - assert top_data[0,0] == 100.0 - ic_package = flopy.mf6.ModflowGwfic(model, strt='initial_heads.txt', - filename='{}.ic'.format(model_name)) - npf_package = ModflowGwfnpf(model, pname='npf_1', save_flows=True, - alternative_cell_averaging='logarithmic', - icelltype=1, k=5.0) + assert top_data[0, 0] == 100.0 + ic_package = flopy.mf6.ModflowGwfic( + model, strt="initial_heads.txt", filename="{}.ic".format(model_name) + ) + npf_package = ModflowGwfnpf( + model, + pname="npf_1", + save_flows=True, + alternative_cell_averaging="logarithmic", + icelltype=1, + k=5.0, + ) # remove package test using .remove_package(name) - assert (model.get_package(npf_package.package_name) is not None) + assert model.get_package(npf_package.package_name) is not None model.remove_package(npf_package.package_name) - assert (model.get_package(npf_package.package_name) is None) + assert model.get_package(npf_package.package_name) is None # remove package test using .remove() - npf_package = ModflowGwfnpf(model, pname='npf_1', save_flows=True, - alternative_cell_averaging='logarithmic', - icelltype=1, k=5.0) + npf_package = ModflowGwfnpf( + model, + pname="npf_1", + save_flows=True, + alternative_cell_averaging="logarithmic", + icelltype=1, + k=5.0, + ) npf_package.remove() - assert (model.get_package(npf_package.package_name) is None) - - npf_package = ModflowGwfnpf(model, save_flows=True, - alternative_cell_averaging='logarithmic', - icelltype=1, k=5.0) - - oc_package = ModflowGwfoc(model, budget_filerecord=[('np001_mod.cbc',)], - head_filerecord=[('np001_mod.hds',)], - saverecord={0: [('HEAD', 'ALL'), - ('BUDGET', 'ALL')], - 1: [('HEAD', 'ALL'), - ('BUDGET', 'ALL')]}, - printrecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) + assert model.get_package(npf_package.package_name) is None + + npf_package = ModflowGwfnpf( + model, + save_flows=True, + alternative_cell_averaging="logarithmic", + icelltype=1, + k=5.0, + ) + + oc_package = ModflowGwfoc( + model, + budget_filerecord=[("np001_mod.cbc",)], + head_filerecord=[("np001_mod.hds",)], + saverecord={ + 0: [("HEAD", "ALL"), ("BUDGET", "ALL")], + 1: [("HEAD", "ALL"), ("BUDGET", "ALL")], + }, + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) oc_package.printrecord.add_transient_key(1) - oc_package.printrecord.set_data([('HEAD', 'ALL'), ('BUDGET', 'ALL')], 1) - - sto_package = ModflowGwfsto(model, save_flows=True, iconvert=1, - ss=0.000001, sy=0.15) + oc_package.printrecord.set_data([("HEAD", "ALL"), ("BUDGET", "ALL")], 1) + + sto_package = ModflowGwfsto( + model, save_flows=True, iconvert=1, ss=0.000001, sy=0.15 + ) + + # verify templates with and without cellid expanded + pkd = ModflowGwfsfr.packagedata.empty(model) + assert pkd.dtype.names[1] == "cellid" + pkd_ex = ModflowGwfsfr.packagedata.empty(model, cellid_expanded=True) + assert pkd_ex.dtype.names[1] == "layer" + assert pkd_ex.dtype.names[2] == "row" + assert pkd_ex.dtype.names[3] == "column" + + pkd_dtype = ModflowGwfsfr.packagedata.dtype(model) + assert pkd_dtype[1][0] == "cellid" + pkd_dtype_ex = ModflowGwfsfr.packagedata.dtype(model, cellid_expanded=True) + assert pkd_dtype_ex[1][0] == "layer" + assert pkd_dtype_ex[2][0] == "row" + assert pkd_dtype_ex[3][0] == "column" # test saving a binary file with list data - well_spd = {0: {'filename': 'wel0.bin', 'binary': True, - 'data': [((0, 0, 4), -2000.0), ((0, 0, 7), -2.0)]}, - 1: None} - wel_package = ModflowGwfwel(model, print_input=True, print_flows=True, - save_flows=True, maxbound=2, - stress_period_data=well_spd) + well_spd = { + 0: { + "filename": "wel0.bin", + "binary": True, + "data": [(0, 0, 4, -2000.0), (0, 0, 7, -2.0)], + }, + 1: None, + } + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=2, + stress_period_data=well_spd, + ) wel_package.stress_period_data.add_transient_key(1) wel_package.stress_period_data.set_data( - {1: {'filename': 'wel.txt', 'factor': 1.0}}) + {1: {"filename": "wel.txt", "factor": 1.0}} + ) # test getting data from a binary file well_data = wel_package.stress_period_data.get_data(0) assert well_data[0][0] == (0, 0, 4) assert well_data[0][1] == -2000.0 - 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=[((0, 0, 0), 80, 'drn_1')]) - drn_package.ts.time_series_namerecord = 'drn_1' - drn_package.ts.interpolation_methodrecord = 'linearend' - - riv_spd = {0: {'filename': 'riv.txt', 'data':[((0, 0, 9), 110, 90.0, - 100.0, 1.0, 2.0, 3.0)]}} - riv_package = ModflowGwfriv(model, print_input=True, print_flows=True, - save_flows=True, maxbound=1, - auxiliary=['var1', 'var2', 'var3'], - stress_period_data=riv_spd) + 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=[((0, 0, 0), 80, "drn_1")], + ) + drn_package.ts.time_series_namerecord = "drn_1" + drn_package.ts.interpolation_methodrecord = "linearend" + + riv_spd = { + 0: { + "filename": "riv.txt", + "data": [((0, 0, 9), 110, 90.0, 100.0, 1.0, 2.0, 3.0)], + } + } + riv_package = ModflowGwfriv( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=1, + auxiliary=["var1", "var2", "var3"], + stress_period_data=riv_spd, + ) riv_data = riv_package.stress_period_data.get_data(0) assert riv_data[0][0] == (0, 0, 9) assert riv_data[0][1] == 110 @@ -230,28 +336,45 @@ def np001(): # verify package look-up pkgs = model.get_package() - assert (len(pkgs) == 9) - pkg = model.get_package('oc') + assert len(pkgs) == 9 + pkg = model.get_package("oc") assert isinstance(pkg, ModflowGwfoc) - pkg = sim.get_package('tdis') + pkg = sim.get_package("tdis") assert isinstance(pkg, ModflowTdis) - pkg = model.get_package('mydispkg') - assert isinstance(pkg, - flopy.mf6.ModflowGwfdis) and \ - pkg.package_name == 'mydispkg' + pkg = model.get_package("mydispkg") + assert ( + isinstance(pkg, flopy.mf6.ModflowGwfdis) + and pkg.package_name == "mydispkg" + ) pkg = model.mydispkg - assert isinstance(pkg, - flopy.mf6.ModflowGwfdis) and \ - pkg.package_name == 'mydispkg' - + assert ( + isinstance(pkg, flopy.mf6.ModflowGwfdis) + and pkg.package_name == "mydispkg" + ) # verify external file contents array_util = PyListUtil() ic_data = ic_package.strt ic_array = ic_data.get_data() - assert array_util.array_comp(ic_array, [[[100.0, 100.0, 100.0, 100.0, - 100.0, 100.0, 100.0, 100.0, - 100.0, 100.0]]]) + assert array_util.array_comp( + ic_array, + [ + [ + [ + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + ] + ] + ], + ) # make folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -266,19 +389,22 @@ def np001(): # get expected results budget_file = os.path.join(os.getcwd(), expected_cbc_file) - budget_obj = bf.CellBudgetFile(budget_file, precision='double') + budget_obj = bf.CellBudgetFile(budget_file, precision="double") budget_frf_valid = np.array( - budget_obj.get_data(text='FLOW-JA-FACE', full3D=True)) + budget_obj.get_data(text="FLOW-JA-FACE", full3D=True) + ) # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'np001_mod.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "np001_mod.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) budget_frf = sim.simulation_data.mfdata[ - (model_name, 'CBC', 'FLOW-JA-FACE')] + (model_name, "CBC", "FLOW-JA-FACE") + ] assert array_util.array_comp(budget_frf_valid, budget_frf) # clean up @@ -286,137 +412,257 @@ def np001(): try: error_occurred = False - well_spd = {0: {'filename': 'wel0.bin', 'binary': True, - 'data': [((0, 0, 4), -2000.0), ((0, 0, 7), -2.0)]}} - wel_package = ModflowGwfwel(model, boundnames=True, - print_input=True, print_flows=True, - save_flows=True, maxbound=2, - stress_period_data=well_spd) + well_spd = { + 0: { + "filename": "wel0.bin", + "binary": True, + "data": [((0, 0, 4), -2000.0), ((0, 0, 7), -2.0)], + } + } + wel_package = ModflowGwfwel( + model, + boundnames=True, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=2, + 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, 'drn_2')]) - npf_package = ModflowGwfnpf(model, save_flows=True, - alternative_cell_averaging='logarithmic', - icelltype=1, k=100001.0, k33=1e-12) + 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, "drn_2"), + ], + ) + 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 + 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 + ) data_invalid = False try: - 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=[((0, 0, 0), 10.0)]) + 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=[((0, 0, 0), 10.0)], + ) except MFDataException: data_invalid = True assert data_invalid + # test -1, -1, -1 cellid + well_spd = {0: [(-1, -1, -1, -2000.0), (0, 0, 7, -2.0)]} + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=2, + stress_period_data=well_spd, + ) + wel_package.write() + mpath = sim.simulation_data.mfpath.get_model_path(model.name) + found_cellid = False + with open(os.path.join(mpath, "np001_mod.wel"), "r") as fd: + for line in fd: + line_lst = line.strip().split() + if ( + len(line) > 2 + and line_lst[0] == "0" + and line_lst[1] == "0" + and line_lst[2] == "0" + ): + found_cellid = True + assert found_cellid + return def np002(): # init paths - test_ex_name = 'np002' - model_name = 'np002_mod' + test_ex_name = "np002" + model_name = "np002_mod" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) - pth_for_mf = os.path.join('..', '..', '..', pth) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) + pth_for_mf = os.path.join("..", "..", "..", pth) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'np002_mod.hds') - expected_cbc_file = os.path.join(expected_output_folder, 'np002_mod.cbc') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "np002_mod.hds") + expected_cbc_file = os.path.join(expected_output_folder, "np002_mod.cbc") # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=run_folder, nocheck=True) + sim = MFSimulation( + sim_name=test_ex_name, + version="mf6", + exe_name=exe_name, + sim_ws=run_folder, + nocheck=True, + ) name = sim.name_file assert name.continue_.get_data() == None assert name.nocheck.get_data() == True assert name.memory_print_option.get_data() == None tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=2, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='ALL', complexity='SIMPLE', - outer_hclose=0.00001, - outer_maximum=50, under_relaxation='NONE', - inner_maximum=30, - inner_hclose=0.00001, linear_acceleration='CG', - preconditioner_levels=7, - preconditioner_drop_tolerance=0.01, - number_orthogonalizations=2) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=2, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="ALL", + complexity="SIMPLE", + outer_hclose=0.00001, + outer_maximum=50, + under_relaxation="NONE", + inner_maximum=30, + inner_hclose=0.00001, + linear_acceleration="CG", + preconditioner_levels=7, + preconditioner_drop_tolerance=0.01, + number_orthogonalizations=2, + ) sim.register_ims_package(ims_package, [model.name]) # get rid of top_data.txt so that a later test does not automatically pass - top_data_file = os.path.join(run_folder, 'top_data.txt') + top_data_file = os.path.join(run_folder, "top_data.txt") if os.path.isfile(top_data_file): os.remove(top_data_file) # test loading data to be stored in a file and loading data from a file # using the "dictionary" input format - top = {'filename': 'top_data.txt', 'factor': 1.0, - 'data': [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, - 100.0, 100.0]} - botm_file = os.path.join(pth_for_mf, 'botm.txt') - botm = {'filename': botm_file, 'factor': 1.0} - dis_package = ModflowGwfdis(model, length_units='FEET', nlay=1, nrow=1, - ncol=10, delr=500.0, delc=500.0, - top=top, botm=botm, idomain=2, - filename='{}.dis'.format(model_name)) - ic_vals = [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, - 100.0] - ic_package = ModflowGwfic(model, strt=ic_vals, - filename='{}.ic'.format(model_name)) - ic_package.strt.store_as_external_file('initial_heads.txt') + top = { + "filename": "top_data.txt", + "factor": 1.0, + "data": [ + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + ], + } + botm_file = os.path.join(pth_for_mf, "botm.txt") + botm = {"filename": botm_file, "factor": 1.0} + dis_package = ModflowGwfdis( + model, + length_units="FEET", + nlay=1, + nrow=1, + ncol=10, + delr=500.0, + delc=500.0, + top=top, + botm=botm, + idomain=2, + filename="{}.dis".format(model_name), + ) + ic_vals = [ + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + 100.0, + ] + ic_package = ModflowGwfic( + model, strt=ic_vals, filename="{}.ic".format(model_name) + ) + ic_package.strt.store_as_external_file("initial_heads.txt") npf_package = ModflowGwfnpf(model, save_flows=True, icelltype=1, k=100.0) - npf_package.k.store_as_external_file('k.bin', binary=True) - oc_package = ModflowGwfoc(model, budget_filerecord=[('np002_mod.cbc',)], - head_filerecord=[('np002_mod.hds',)], - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) + npf_package.k.store_as_external_file("k.bin", binary=True) + oc_package = ModflowGwfoc( + model, + budget_filerecord=[("np002_mod.cbc",)], + head_filerecord=[("np002_mod.hds",)], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) oc_package.saverecord.add_transient_key(1) - oc_package.saverecord.set_data([('HEAD', 'ALL'), ('BUDGET', 'ALL')], 1) + oc_package.saverecord.set_data([("HEAD", "ALL"), ("BUDGET", "ALL")], 1) oc_package.printrecord.add_transient_key(1) - oc_package.printrecord.set_data([('HEAD', 'ALL'), ('BUDGET', 'ALL')], 1) - - sto_package = ModflowGwfsto(model, save_flows=True, iconvert=1, - ss=0.000001, sy=0.15) - - hfb_package = ModflowGwfhfb(model, print_input=True, maxhfb=1, - stress_period_data=[((0, 0, 3), (0, 0, 4), - 0.00001)]) - chd_package = ModflowGwfchd(model, print_input=True, print_flows=True, - maxbound=1, stress_period_data=[((0, 0, 0), - 65.0)]) - ghb_package = ModflowGwfghb(model, print_input=True, print_flows=True, - maxbound=1, stress_period_data=[((0, 0, 9), - 125.0, 60.0)]) - rch_package = ModflowGwfrch(model, print_input=True, print_flows=True, - maxbound=2, - stress_period_data=[((0, 0, 3), 0.02), - ((0, 0, 6), 0.1)]) + oc_package.printrecord.set_data([("HEAD", "ALL"), ("BUDGET", "ALL")], 1) + + sto_package = ModflowGwfsto( + model, save_flows=True, iconvert=1, ss=0.000001, sy=0.15 + ) + + hfb_package = ModflowGwfhfb( + model, + print_input=True, + maxhfb=1, + stress_period_data=[((0, 0, 3), (0, 0, 4), 0.00001)], + ) + chd_package = ModflowGwfchd( + model, + print_input=True, + print_flows=True, + maxbound=1, + stress_period_data=[((0, 0, 0), 65.0)], + ) + ghb_package = ModflowGwfghb( + model, + print_input=True, + print_flows=True, + maxbound=1, + stress_period_data=[((0, 0, 9), 125.0, 60.0)], + ) + rch_package = ModflowGwfrch( + model, + print_input=True, + print_flows=True, + maxbound=2, + stress_period_data=[((0, 0, 3), 0.02), ((0, 0, 6), 0.1)], + ) # write simulation to new location sim.write_simulation() - assert(os.path.isfile(top_data_file)) + assert os.path.isfile(top_data_file) if run: # run simulation @@ -424,30 +670,33 @@ def np002(): sim2 = MFSimulation.load(sim_ws=run_folder) model_ = sim2.get_model(model_name) - npf_package = model_.get_package('npf') + npf_package = model_.get_package("npf") k = npf_package.k.array # get expected results budget_file = os.path.join(os.getcwd(), expected_cbc_file) - budget_obj = bf.CellBudgetFile(budget_file, precision='double') + budget_obj = bf.CellBudgetFile(budget_file, precision="double") budget_frf_valid = np.array( - budget_obj.get_data(text='FLOW JA FACE ', full3D=True)) + budget_obj.get_data(text="FLOW JA FACE ", full3D=True) + ) # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'np002_mod.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, - files2=head_new, outfile=outfile) + head_new = os.path.join(run_folder, "np002_mod.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) array_util = PyListUtil() budget_frf = sim.simulation_data.mfdata[ - (model_name, 'CBC', 'FLOW-JA-FACE')] + (model_name, "CBC", "FLOW-JA-FACE") + ] assert array_util.array_comp(budget_frf_valid, budget_frf) # verify external text file was written correctly - ext_file_path = os.path.join(run_folder, 'initial_heads.txt') - fd = open(ext_file_path, 'r') + ext_file_path = os.path.join(run_folder, "initial_heads.txt") + fd = open(ext_file_path, "r") line = fd.readline() line_array = line.split() assert len(ic_vals) == len(line_array) @@ -459,116 +708,188 @@ def np002(): 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)]) + 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 + 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 def test021_twri(): # init paths - test_ex_name = 'test021_twri' - model_name = 'twri' + test_ex_name = "test021_twri" + model_name = "twri" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'twri.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "twri.hds") # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=test_ex_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) tdis_rc = [(86400.0, 1, 1.0)] - tdis_package = ModflowTdis(sim, time_units='SECONDS', nper=1, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='SUMMARY', outer_hclose=0.0001, - outer_maximum=500, under_relaxation='NONE', - inner_maximum=100, - inner_hclose=0.0001, rcloserecord=0.001, - linear_acceleration='CG', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.97) + tdis_package = ModflowTdis( + sim, time_units="SECONDS", nper=1, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + outer_hclose=0.0001, + outer_maximum=500, + under_relaxation="NONE", + inner_maximum=100, + inner_hclose=0.0001, + rcloserecord=0.001, + linear_acceleration="CG", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.97, + ) sim.register_ims_package(ims_package, [model.name]) # build top binary data - text = 'TOP' - fname = 'top.bin' + text = "TOP" + fname = "top.bin" nrow = 15 ncol = 15 pth = os.path.join(sim.simulation_data.mfpath.get_sim_path(), fname) - f = open(pth, 'wb') - header = flopy.utils.BinaryHeader.create(bintype='HEAD', - precision='double', - text=text, - nrow=nrow, - ncol=ncol, - ilay=1, pertim=1.0, - totim=1.0, kstp=1, - kper=1) - flopy.utils.Util2d.write_bin((nrow, ncol), f, - np.full((nrow, ncol), 200.0, - dtype=np.float64), header_data=header) + f = open(pth, "wb") + header = flopy.utils.BinaryHeader.create( + bintype="HEAD", + precision="double", + text=text, + nrow=nrow, + ncol=ncol, + ilay=1, + pertim=1.0, + totim=1.0, + kstp=1, + kper=1, + ) + flopy.utils.Util2d.write_bin( + (nrow, ncol), + f, + np.full((nrow, ncol), 200.0, dtype=np.float64), + header_data=header, + ) f.close() - top = {'factor': 1., 'filename': fname, 'data': None, 'binary': True, - 'iprn': 1} - - dis_package = flopy.mf6.ModflowGwfdis(model, nlay=3, nrow=15, ncol=15, - delr=5000.0, delc=5000.0, - top=top, botm=[-200, -300, -450], - filename='{}.dis'.format(model_name)) - strt = [{'filename': 'strt.txt', 'factor': 1.0, 'data': 0.0}, - {'filename': 'strt2.bin', 'factor': 1.0, 'data': 1.0, - 'binary': 'True'}, 2.0] - ic_package = ModflowGwfic(model, strt=strt, - filename='{}.ic'.format(model_name)) - npf_package = ModflowGwfnpf(model, save_flows=True, perched=True, - cvoptions='dewatered', - icelltype=[1, 0, 0], k=[0.001, 0.0001, 0.0002], - k33=0.00000002) - oc_package = ModflowGwfoc(model, budget_filerecord='twri.cbc', - head_filerecord='twri.hds', - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL')]) + top = { + "factor": 1.0, + "filename": fname, + "data": None, + "binary": True, + "iprn": 1, + } + + dis_package = flopy.mf6.ModflowGwfdis( + model, + nlay=3, + nrow=15, + ncol=15, + delr=5000.0, + delc=5000.0, + top=top, + botm=[-200, -300, -450], + filename="{}.dis".format(model_name), + ) + strt = [ + {"filename": "strt.txt", "factor": 1.0, "data": 0.0}, + { + "filename": "strt2.bin", + "factor": 1.0, + "data": 1.0, + "binary": "True", + }, + 2.0, + ] + ic_package = ModflowGwfic( + model, strt=strt, filename="{}.ic".format(model_name) + ) + npf_package = ModflowGwfnpf( + model, + save_flows=True, + perched=True, + cvoptions="dewatered", + icelltype=[1, 0, 0], + k=[0.001, 0.0001, 0.0002], + k33=0.00000002, + ) + oc_package = ModflowGwfoc( + model, + budget_filerecord="twri.cbc", + head_filerecord="twri.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL")], + ) # build stress_period_data for chd package stress_period_data = [] for layer in range(0, 2): for row in range(0, 15): stress_period_data.append(((layer, row, 0), 0.0)) - chd_package = ModflowGwfchd(model, print_input=True, print_flows=True, - save_flows=True, maxbound=100, - stress_period_data=stress_period_data) + chd_package = ModflowGwfchd( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=100, + stress_period_data=stress_period_data, + ) # build stress_period_data for drn package - conc = np.ones((15, 15), dtype=np.float) * 35. + conc = np.ones((15, 15), dtype=np.float) * 35.0 auxdata = {0: [6, conc]} stress_period_data = [] drn_heads = [0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0] for col, head in zip(range(1, 10), drn_heads): - stress_period_data.append(((0, 7, col), head, 1.0, - 'name_{}'.format(col))) - drn_package = ModflowGwfdrn(model, print_input=True, print_flows=True, - save_flows=True, maxbound=9, boundnames=True, - stress_period_data=stress_period_data) - rch_package = ModflowGwfrcha(model, readasarrays=True, fixed_cell=True, - recharge={0: 0.00000003}, - auxiliary=[('iface', 'conc')], aux=auxdata) + stress_period_data.append( + ((0, 7, col), head, 1.0, "name_{}".format(col)) + ) + drn_package = ModflowGwfdrn( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=9, + boundnames=True, + stress_period_data=stress_period_data, + ) + rch_package = ModflowGwfrcha( + model, + readasarrays=True, + fixed_cell=True, + recharge={0: 0.00000003}, + auxiliary=[("iface", "conc")], + aux=auxdata, + ) aux = rch_package.aux.get_data() @@ -578,9 +899,14 @@ def test021_twri(): cols = [10, 5, 11, 7, 9, 11, 13, 7, 9, 11, 13, 7, 9, 11, 13] for layer, row, col in zip(layers, rows, cols): stress_period_data.append(((layer, row, col), -5.0)) - wel_package = ModflowGwfwel(model, print_input=True, print_flows=True, - save_flows=True, maxbound=15, - stress_period_data=stress_period_data) + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=15, + stress_period_data=stress_period_data, + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -593,21 +919,22 @@ def test021_twri(): sim2 = MFSimulation.load(sim_ws=run_folder) model2 = sim2.get_model() - ic2 = model2.get_package('ic') + ic2 = model2.get_package("ic") strt2 = ic2.strt.get_data() - drn2 = model2.get_package('drn') + drn2 = model2.get_package("drn") drn_spd = drn2.stress_period_data.get_data() - assert(strt2[0,0,0] == 0.0) - assert(strt2[1,0,0] == 1.0) - assert(strt2[2,0,0] == 2.0) - assert(drn_spd[0][1][3] == 'name_2') + assert strt2[0, 0, 0] == 0.0 + assert strt2[1, 0, 0] == 1.0 + assert strt2[2, 0, 0] == 2.0 + assert drn_spd[0][1][3] == "name_2" # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'twri.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "twri.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # clean up sim.delete_output_files() @@ -617,102 +944,142 @@ def test021_twri(): def test005_advgw_tidal(): # init paths - test_ex_name = 'test005_advgw_tidal' - model_name = 'AdvGW_tidal' + test_ex_name = "test005_advgw_tidal" + model_name = "AdvGW_tidal" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, - 'AdvGW_tidal.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join( + expected_output_folder, "AdvGW_tidal.hds" + ) # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=test_ex_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) # test tdis package deletion - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=1, - perioddata=[(2.0, 2, 1.0)]) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=1, perioddata=[(2.0, 2, 1.0)] + ) sim.remove_package(tdis_package.package_type) - tdis_rc = [(1.0, 1, 1.0), (10.0, 120, 1.0), (10.0, 120, 1.0), - (10.0, 120, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=4, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='SUMMARY', complexity='SIMPLE', - outer_hclose=0.0001, - outer_maximum=500, under_relaxation='NONE', - inner_maximum=100, - inner_hclose=0.0001, rcloserecord=0.001, - linear_acceleration='CG', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.97) + tdis_rc = [ + (1.0, 1, 1.0), + (10.0, 120, 1.0), + (10.0, 120, 1.0), + (10.0, 120, 1.0), + ] + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=4, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + complexity="SIMPLE", + outer_hclose=0.0001, + outer_maximum=500, + under_relaxation="NONE", + inner_maximum=100, + inner_hclose=0.0001, + rcloserecord=0.001, + linear_acceleration="CG", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.97, + ) sim.register_ims_package(ims_package, [model.name]) bot_data = [-100 for x in range(150)] - dis_package = ModflowGwfdis(model, nlay=3, nrow=15, ncol=10, delr=500.0, - delc=500.0, - top=50.0, botm=[5.0, -10.0, {'factor': 1.0, - 'data': bot_data}], - filename='{}.dis'.format(model_name)) - ic_package = ModflowGwfic(model, strt=50.0, - filename='{}.ic'.format(model_name)) - npf_package = ModflowGwfnpf(model, save_flows=True, icelltype=[1, 0, 0], - k=[5.0, 0.1, 4.0], - k33=[0.5, 0.005, 0.1]) - oc_package = ModflowGwfoc(model, budget_filerecord='AdvGW_tidal.cbc', - head_filerecord='AdvGW_tidal.hds', - headprintrecord=[('COLUMNS', 10, 'WIDTH', 15, - 'DIGITS', 6, 'GENERAL')], - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'FIRST'), ('HEAD', 'LAST'), - ('BUDGET', 'LAST')]) + dis_package = ModflowGwfdis( + model, + nlay=3, + nrow=15, + ncol=10, + delr=500.0, + delc=500.0, + top=50.0, + botm=[5.0, -10.0, {"factor": 1.0, "data": bot_data}], + filename="{}.dis".format(model_name), + ) + ic_package = ModflowGwfic( + model, strt=50.0, filename="{}.ic".format(model_name) + ) + npf_package = ModflowGwfnpf( + model, + save_flows=True, + icelltype=[1, 0, 0], + k=[5.0, 0.1, 4.0], + k33=[0.5, 0.005, 0.1], + ) + oc_package = ModflowGwfoc( + model, + budget_filerecord="AdvGW_tidal.cbc", + head_filerecord="AdvGW_tidal.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "FIRST"), ("HEAD", "LAST"), ("BUDGET", "LAST")], + ) # test empty sy_template = ModflowGwfsto.sy.empty(model, True) for layer in range(0, 3): - sy_template[layer]['data'] = 0.2 - layer_storage_types = [DataStorageType.internal_array, - DataStorageType.internal_constant, - DataStorageType.internal_array] - ss_template = ModflowGwfsto.ss.empty(model, True, layer_storage_types, - 0.000001) - sto_package = ModflowGwfsto(model, save_flows=True, iconvert=1, - ss=ss_template, sy=sy_template, - steady_state={0: True}, - transient={1: True}) + sy_template[layer]["data"] = 0.2 + layer_storage_types = [ + DataStorageType.internal_array, + DataStorageType.internal_constant, + DataStorageType.internal_array, + ] + ss_template = ModflowGwfsto.ss.empty( + model, True, layer_storage_types, 0.000001 + ) + sto_package = ModflowGwfsto( + model, + save_flows=True, + iconvert=1, + ss=ss_template, + sy=sy_template, + steady_state={0: True}, + transient={1: True}, + ) # wel, evt, ghb, obs, riv, rch, ts # well package # test empty with aux vars, bound names, and time series - period_two = ModflowGwfwel.stress_period_data.empty(model, maxbound=3, - aux_vars=['var1', - 'var2', - 'var3'], - boundnames=True, - timeseries=True) + period_two = ModflowGwfwel.stress_period_data.empty( + model, + maxbound=3, + aux_vars=["var1", "var2", "var3"], + boundnames=True, + timeseries=True, + ) period_two[0][0] = ((0, 11, 2), -50.0, -1, -2, -3, None) - period_two[0][1] = ((2, 4, 7), 'well_1_rate', 1, 2, 3, 'well_1') - period_two[0][2] = ((2, 3, 2), 'well_2_rate', 4, 5, 6, 'well_2') - period_three = ModflowGwfwel.stress_period_data.empty(model, maxbound=2, - aux_vars=['var1', - 'var2', - 'var3'], - boundnames=True, - timeseries=True) - period_three[0][0] = ((2, 3, 2), 'well_2_rate', 1, 2, 3, 'well_2') - period_three[0][1] = ((2, 4, 7), 'well_1_rate', 4, 5, 6, 'well_1') - period_four = ModflowGwfwel.stress_period_data.empty(model, maxbound=5, - aux_vars=['var1', - 'var2', - 'var3'], - boundnames=True, - timeseries=True) - period_four[0][0] = ((2, 4, 7), 'well_1_rate', 1, 2, 3, 'well_1') - period_four[0][1] = ((2, 3, 2), 'well_2_rate', 4, 5, 6, 'well_2') + period_two[0][1] = ((2, 4, 7), "well_1_rate", 1, 2, 3, "well_1") + period_two[0][2] = ((2, 3, 2), "well_2_rate", 4, 5, 6, "well_2") + period_three = ModflowGwfwel.stress_period_data.empty( + model, + maxbound=2, + aux_vars=["var1", "var2", "var3"], + boundnames=True, + timeseries=True, + ) + period_three[0][0] = ((2, 3, 2), "well_2_rate", 1, 2, 3, "well_2") + period_three[0][1] = ((2, 4, 7), "well_1_rate", 4, 5, 6, "well_1") + period_four = ModflowGwfwel.stress_period_data.empty( + model, + maxbound=5, + aux_vars=["var1", "var2", "var3"], + boundnames=True, + timeseries=True, + ) + period_four[0][0] = ((2, 4, 7), "well_1_rate", 1, 2, 3, "well_1") + period_four[0][1] = ((2, 3, 2), "well_2_rate", 4, 5, 6, "well_2") period_four[0][2] = ((0, 11, 2), -10.0, 7, 8, 9, None) period_four[0][3] = ((0, 2, 4), -20.0, 17, 18, 19, None) period_four[0][4] = ((0, 13, 5), -40.0, 27, 28, 29, None) @@ -721,40 +1088,70 @@ def test005_advgw_tidal(): stress_period_data[2] = period_three[0] stress_period_data[3] = period_four[0] # well ts package - timeseries = [(0.0, 0.0, 0.0, 0.0), - (1.0, -200.0, 0.0, -100.0), - (11.0, -1800.0, -500.0, -200.0), - (21.0, -200.0, -400.0, -300.0), - (31.0, 0.0, -600.0, -400.0)] - ts_dict = {'filename': 'well-rates.ts', 'timeseries': timeseries, - 'time_series_namerecord': [('well_1_rate', 'well_2_rate', - 'well_3_rate')], - 'interpolation_methodrecord': [('stepwise', 'stepwise', - 'stepwise')]} + timeseries = [ + (0.0, 0.0, 0.0, 0.0), + (1.0, -200.0, 0.0, -100.0), + (11.0, -1800.0, -500.0, -200.0), + (21.0, -200.0, -400.0, -300.0), + (31.0, 0.0, -600.0, -400.0), + ] + ts_dict = { + "filename": "well-rates.ts", + "timeseries": timeseries, + "time_series_namerecord": [ + ("well_1_rate", "well_2_rate", "well_3_rate") + ], + "interpolation_methodrecord": [("stepwise", "stepwise", "stepwise")], + } # test removing package with child packages - wel_package = ModflowGwfwel(model, print_input=True, print_flows=True, - auxiliary=[('var1', 'var2', 'var3')], - maxbound=5, - stress_period_data=stress_period_data, - boundnames=True, save_flows=True, - timeseries=ts_dict) + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + auxiliary=[("var1", "var2", "var3")], + maxbound=5, + stress_period_data=stress_period_data, + boundnames=True, + save_flows=True, + timeseries=ts_dict, + ) wel_package.remove() - wel_package = ModflowGwfwel(model, print_input=True, print_flows=True, - auxiliary=[('var1', 'var2', 'var3')], - maxbound=5, - stress_period_data=stress_period_data, - boundnames=True, save_flows=True, - timeseries=ts_dict) + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + auxiliary=[("var1", "var2", "var3")], + maxbound=5, + stress_period_data=stress_period_data, + boundnames=True, + save_flows=True, + timeseries=ts_dict, + ) # test empty evt_period = ModflowGwfevt.stress_period_data.empty(model, 150, nseg=3) for col in range(0, 10): for row in range(0, 15): evt_period[0][col * 15 + row] = ( - ((0, row, col), 50.0, 0.0004, 10.0, 0.2, 0.5, 0.3, 0.1, None)) - evt_package = ModflowGwfevt(model, print_input=True, print_flows=True, - save_flows=True, maxbound=150, - nseg=3, stress_period_data=evt_period) + (0, row, col), + 50.0, + 0.0004, + 10.0, + 0.2, + 0.5, + 0.3, + 0.1, + None, + ) + evt_package = ModflowGwfevt( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=150, + nseg=3, + stress_period_data=evt_period, + ) # build ghb ghb_period = {} @@ -762,94 +1159,133 @@ def test005_advgw_tidal(): for layer, cond in zip(range(1, 3), [15.0, 1500.0]): for row in range(0, 15): ghb_period_array.append( - ((layer, row, 9), 'tides', cond, 'Estuary-L2')) + ((layer, row, 9), "tides", cond, "Estuary-L2") + ) ghb_period[0] = ghb_period_array # build ts ghb ts_recarray = [] - fd = open(os.path.join(pth, 'tides.txt'), 'r') + fd = open(os.path.join(pth, "tides.txt"), "r") for line in fd: - line_list = line.strip().split(',') + line_list = line.strip().split(",") ts_recarray.append((float(line_list[0]), float(line_list[1]))) - ts_package_dict = {'filename':'tides.ts', - 'timeseries':ts_recarray, - 'time_series_namerecord':'tides', - 'interpolation_methodrecord':'linear'} - - obs_dict = {('ghb_obs.csv', 'binary'): [('ghb-2-6-10', 'GHB', (1, 5, 9)), - ('ghb-3-6-10', 'GHB', (2, 5, 9))], - 'ghb_flows.csv': [('Estuary2', 'GHB', 'Estuary-L2'), - ('Estuary3', 'GHB', 'Estuary-L3')], - 'filename': 'AdvGW_tidal.ghb.obs', 'digits': 10, - 'print_input': True} - - ghb_package = ModflowGwfghb(model, print_input=True, print_flows=True, - save_flows=True, boundnames=True, - timeseries=ts_package_dict, - observations=obs_dict, - maxbound=30, stress_period_data=ghb_period) + ts_package_dict = { + "filename": "tides.ts", + "timeseries": ts_recarray, + "time_series_namerecord": "tides", + "interpolation_methodrecord": "linear", + } + + obs_dict = { + ("ghb_obs.csv", "binary"): [ + ("ghb-2-6-10", "GHB", (1, 5, 9)), + ("ghb-3-6-10", "GHB", (2, 5, 9)), + ], + "ghb_flows.csv": [ + ("Estuary2", "GHB", "Estuary-L2"), + ("Estuary3", "GHB", "Estuary-L3"), + ], + "filename": "AdvGW_tidal.ghb.obs", + "digits": 10, + "print_input": True, + } + + ghb_package = ModflowGwfghb( + model, + print_input=True, + print_flows=True, + save_flows=True, + boundnames=True, + timeseries=ts_package_dict, + observations=obs_dict, + maxbound=30, + stress_period_data=ghb_period, + ) riv_period = {} - riv_period_array = [((0, 2, 0), 'river_stage_1', 1001.0, 35.9, None), - ((0, 3, 1), 'river_stage_1', 1002.0, 35.8, None), - ((0, 4, 2), 'river_stage_1', 1003.0, 35.7, None), - ((0, 4, 3), 'river_stage_1', 1004.0, 35.6, None), - ((0, 5, 4), 'river_stage_1', 1005.0, 35.5, None), - ((0, 5, 5), 'river_stage_1', 1006.0, 35.4, 'riv1_c6'), - ((0, 5, 6), 'river_stage_1', 1007.0, 35.3, 'riv1_c7'), - ((0, 4, 7), 'river_stage_1', 1008.0, 35.2, None), - ((0, 4, 8), 'river_stage_1', 1009.0, 35.1, None), - ((0, 4, 9), 'river_stage_1', 1010.0, 35.0, None), - ((0, 9, 0), 'river_stage_2', 1001.0, 36.9, - 'riv2_upper'), - ((0, 8, 1), 'river_stage_2', 1002.0, 36.8, - 'riv2_upper'), - ((0, 7, 2), 'river_stage_2', 1003.0, 36.7, - 'riv2_upper'), - ((0, 6, 3), 'river_stage_2', 1004.0, 36.6, None), - ((0, 6, 4), 'river_stage_2', 1005.0, 36.5, None), - ((0, 5, 5), 'river_stage_2', 1006.0, 36.4, 'riv2_c6'), - ((0, 5, 6), 'river_stage_2', 1007.0, 36.3, 'riv2_c7'), - ((0, 6, 7), 'river_stage_2', 1008.0, 36.2, None), - ((0, 6, 8), 'river_stage_2', 1009.0, 36.1), - ((0, 6, 9), 'river_stage_2', 1010.0, 36.0)] + riv_period_array = [ + ((0, 2, 0), "river_stage_1", 1001.0, 35.9, None), + ((0, 3, 1), "river_stage_1", 1002.0, 35.8, None), + ((0, 4, 2), "river_stage_1", 1003.0, 35.7, None), + ((0, 4, 3), "river_stage_1", 1004.0, 35.6, None), + ((0, 5, 4), "river_stage_1", 1005.0, 35.5, None), + ((0, 5, 5), "river_stage_1", 1006.0, 35.4, "riv1_c6"), + ((0, 5, 6), "river_stage_1", 1007.0, 35.3, "riv1_c7"), + ((0, 4, 7), "river_stage_1", 1008.0, 35.2, None), + ((0, 4, 8), "river_stage_1", 1009.0, 35.1, None), + ((0, 4, 9), "river_stage_1", 1010.0, 35.0, None), + ((0, 9, 0), "river_stage_2", 1001.0, 36.9, "riv2_upper"), + ((0, 8, 1), "river_stage_2", 1002.0, 36.8, "riv2_upper"), + ((0, 7, 2), "river_stage_2", 1003.0, 36.7, "riv2_upper"), + ((0, 6, 3), "river_stage_2", 1004.0, 36.6, None), + ((0, 6, 4), "river_stage_2", 1005.0, 36.5, None), + ((0, 5, 5), "river_stage_2", 1006.0, 36.4, "riv2_c6"), + ((0, 5, 6), "river_stage_2", 1007.0, 36.3, "riv2_c7"), + ((0, 6, 7), "river_stage_2", 1008.0, 36.2, None), + ((0, 6, 8), "river_stage_2", 1009.0, 36.1), + ((0, 6, 9), "river_stage_2", 1010.0, 36.0), + ] riv_period[0] = riv_period_array # riv time series - ts_data = [(0.0, 40.0, 41.0), (1.0, 41.0, 41.5), (2.0, 43.0, 42.0), - (3.0, 45.0, 42.8), (4.0, 44.0, 43.0), - (6.0, 43.0, 43.1), (9.0, 42.0, 42.4), (11.0, 41.0, 41.5), - (31.0, 40.0, 41.0)] - ts_dict = {'filename': 'river_stages.ts', 'timeseries': ts_data, - 'time_series_namerecord': [('river_stage_1', 'river_stage_2')], - 'interpolation_methodrecord': [('linear', 'stepwise')]} + ts_data = [ + (0.0, 40.0, 41.0), + (1.0, 41.0, 41.5), + (2.0, 43.0, 42.0), + (3.0, 45.0, 42.8), + (4.0, 44.0, 43.0), + (6.0, 43.0, 43.1), + (9.0, 42.0, 42.4), + (11.0, 41.0, 41.5), + (31.0, 40.0, 41.0), + ] + ts_dict = { + "filename": "river_stages.ts", + "timeseries": ts_data, + "time_series_namerecord": [("river_stage_1", "river_stage_2")], + "interpolation_methodrecord": [("linear", "stepwise")], + } # riv obs - obs_dict = {'riv_obs.csv': [('rv1-3-1', 'RIV', (0, 2, 0)), - ('rv1-4-2', 'RIV', (0, 3, 1)), - ('rv1-5-3', 'RIV', (0, 4, 2)), - ('rv1-5-4', 'RIV', (0, 4, 3)), - ('rv1-6-5', 'RIV', (0, 5, 4)), - ('rv1-c6', 'RIV', 'riv1_c6'), - ('rv1-c7', 'RIV', 'riv1_c7'), - ('rv2-upper', 'RIV', 'riv2_upper'), - ('rv-2-7-4', 'RIV', (0, 6, 3)), - ('rv2-8-5', 'RIV', (0, 6, 4)), - ('rv-2-9-6', 'RIV', (0, 5, 5,))], - 'riv_flowsA.csv': [('riv1-3-1', 'RIV', (0, 2, 0)), - ('riv1-4-2', 'RIV', (0, 3, 1)), - ('riv1-5-3', 'RIV', (0, 4, 2))], - 'riv_flowsB.csv': [('riv2-10-1', 'RIV', (0, 9, 0)), - ('riv-2-9-2', 'RIV', (0, 8, 1)), - ('riv2-8-3', 'RIV', (0, 7, 2))], - 'filename': 'AdvGW_tidal.riv.obs', 'digits': 10, - 'print_input': True} - - riv_package = ModflowGwfriv(model, print_input=True, print_flows=True, - save_flows=True, - boundnames=True, - timeseries=ts_dict, - maxbound=20, stress_period_data=riv_period, - observations=obs_dict) + obs_dict = { + "riv_obs.csv": [ + ("rv1-3-1", "RIV", (0, 2, 0)), + ("rv1-4-2", "RIV", (0, 3, 1)), + ("rv1-5-3", "RIV", (0, 4, 2)), + ("rv1-5-4", "RIV", (0, 4, 3)), + ("rv1-6-5", "RIV", (0, 5, 4)), + ("rv1-c6", "RIV", "riv1_c6"), + ("rv1-c7", "RIV", "riv1_c7"), + ("rv2-upper", "RIV", "riv2_upper"), + ("rv-2-7-4", "RIV", (0, 6, 3)), + ("rv2-8-5", "RIV", (0, 6, 4)), + ("rv-2-9-6", "RIV", (0, 5, 5,)), + ], + "riv_flowsA.csv": [ + ("riv1-3-1", "RIV", (0, 2, 0)), + ("riv1-4-2", "RIV", (0, 3, 1)), + ("riv1-5-3", "RIV", (0, 4, 2)), + ], + "riv_flowsB.csv": [ + ("riv2-10-1", "RIV", (0, 9, 0)), + ("riv-2-9-2", "RIV", (0, 8, 1)), + ("riv2-8-3", "RIV", (0, 7, 2)), + ], + "filename": "AdvGW_tidal.riv.obs", + "digits": 10, + "print_input": True, + } + + riv_package = ModflowGwfriv( + model, + print_input=True, + print_flows=True, + save_flows=True, + boundnames=True, + timeseries=ts_dict, + maxbound=20, + stress_period_data=riv_period, + observations=obs_dict, + ) rch1_period = {} rch1_period_array = [] @@ -860,63 +1296,103 @@ def test005_advgw_tidal(): else: col_max = 6 for col in range(0, col_max): - if (row == 3 and col == 5) or (row == 2 and col == 4) or ( - row == 1 and col == 3) or (row == 0 and col == 2): + if ( + (row == 3 and col == 5) + or (row == 2 and col == 4) + or (row == 1 and col == 3) + or (row == 0 and col == 2) + ): mult = 0.5 else: mult = 1.0 if row == 0 and col == 0: - bnd = 'rch-1-1' + bnd = "rch-1-1" elif row == 0 and col == 1: - bnd = 'rch-1-2' + bnd = "rch-1-2" elif row == 1 and col == 2: - bnd = 'rch-2-3' + bnd = "rch-2-3" else: bnd = None - rch1_period_array.append(((0, row, col), 'rch_1', mult, bnd)) + rch1_period_array.append(((0, row, col), "rch_1", mult, bnd)) rch1_period[0] = rch1_period_array - rch1_package = ModflowGwfrch(model, filename='AdvGW_tidal_1.rch', - pname='rch_1', fixed_cell=True, - auxiliary='MULTIPLIER', - auxmultname='MULTIPLIER', - print_input=True, print_flows=True, - save_flows=True, boundnames=True, - maxbound=84, stress_period_data=rch1_period) - ts_data = [(0.0, 0.0015), (1.0, 0.0010), (11.0, 0.0015), - (21.0, 0.0025), (31.0, 0.0015)] - rch1_package.ts.initialize(timeseries=ts_data, - filename='recharge_rates_1.ts', - time_series_namerecord='rch_1', - interpolation_methodrecord='stepwise') + rch1_package = ModflowGwfrch( + model, + filename="AdvGW_tidal_1.rch", + pname="rch_1", + fixed_cell=True, + auxiliary="MULTIPLIER", + auxmultname="MULTIPLIER", + print_input=True, + print_flows=True, + save_flows=True, + boundnames=True, + maxbound=84, + stress_period_data=rch1_period, + ) + ts_data = [ + (0.0, 0.0015), + (1.0, 0.0010), + (11.0, 0.0015), + (21.0, 0.0025), + (31.0, 0.0015), + ] + rch1_package.ts.initialize( + timeseries=ts_data, + filename="recharge_rates_1.ts", + time_series_namerecord="rch_1", + interpolation_methodrecord="stepwise", + ) rch2_period = {} - rch2_period_array = [((0, 0, 2), 'rch_2', 0.5), ((0, 0, 3), 'rch_2', 1.0), - ((0, 0, 4), 'rch_2', 1.0), - ((0, 0, 5), 'rch_2', 1.0), ((0, 0, 6), 'rch_2', 1.0), - ((0, 0, 7), 'rch_2', 1.0), - ((0, 0, 8), 'rch_2', 1.0), ((0, 0, 9), 'rch_2', 0.5), - ((0, 1, 3), 'rch_2', 0.5), - ((0, 1, 4), 'rch_2', 1.0), ((0, 1, 5), 'rch_2', 1.0), - ((0, 1, 6), 'rch_2', 1.0), - ((0, 1, 7), 'rch_2', 1.0), ((0, 1, 8), 'rch_2', 0.5), - ((0, 2, 4), 'rch_2', 0.5), - ((0, 2, 5), 'rch_2', 1.0), ((0, 2, 6), 'rch_2', 1.0), - ((0, 2, 7), 'rch_2', 0.5), - ((0, 3, 5), 'rch_2', 0.5), ((0, 3, 6), 'rch_2', 0.5)] + rch2_period_array = [ + ((0, 0, 2), "rch_2", 0.5), + ((0, 0, 3), "rch_2", 1.0), + ((0, 0, 4), "rch_2", 1.0), + ((0, 0, 5), "rch_2", 1.0), + ((0, 0, 6), "rch_2", 1.0), + ((0, 0, 7), "rch_2", 1.0), + ((0, 0, 8), "rch_2", 1.0), + ((0, 0, 9), "rch_2", 0.5), + ((0, 1, 3), "rch_2", 0.5), + ((0, 1, 4), "rch_2", 1.0), + ((0, 1, 5), "rch_2", 1.0), + ((0, 1, 6), "rch_2", 1.0), + ((0, 1, 7), "rch_2", 1.0), + ((0, 1, 8), "rch_2", 0.5), + ((0, 2, 4), "rch_2", 0.5), + ((0, 2, 5), "rch_2", 1.0), + ((0, 2, 6), "rch_2", 1.0), + ((0, 2, 7), "rch_2", 0.5), + ((0, 3, 5), "rch_2", 0.5), + ((0, 3, 6), "rch_2", 0.5), + ] rch2_period[0] = rch2_period_array - rch2_package = ModflowGwfrch(model, filename='AdvGW_tidal_2.rch', - pname='rch_2', fixed_cell=True, - auxiliary='MULTIPLIER', - auxmultname='MULTIPLIER', - print_input=True, print_flows=True, - save_flows=True, - maxbound=20, stress_period_data=rch2_period) - ts_data = [(0.0, 0.0016), (1.0, 0.0018), (11.0, 0.0019), - (21.0, 0.0016), (31.0, 0.0018)] - rch2_package.ts.initialize(timeseries=ts_data, - filename='recharge_rates_2.ts', - time_series_namerecord='rch_2', - interpolation_methodrecord='linear') + rch2_package = ModflowGwfrch( + model, + filename="AdvGW_tidal_2.rch", + pname="rch_2", + fixed_cell=True, + auxiliary="MULTIPLIER", + auxmultname="MULTIPLIER", + print_input=True, + print_flows=True, + save_flows=True, + maxbound=20, + stress_period_data=rch2_period, + ) + ts_data = [ + (0.0, 0.0016), + (1.0, 0.0018), + (11.0, 0.0019), + (21.0, 0.0016), + (31.0, 0.0018), + ] + rch2_package.ts.initialize( + timeseries=ts_data, + filename="recharge_rates_2.ts", + time_series_namerecord="rch_2", + interpolation_methodrecord="linear", + ) rch3_period = {} rch3_period_array = [] @@ -927,27 +1403,43 @@ def test005_advgw_tidal(): else: col_min = 6 for col in range(col_min, 10): - if (row == 0 and col == 9) or (row == 1 and col == 8) or ( - row == 2 and col == 7) or (row == 3 and col == 6): + if ( + (row == 0 and col == 9) + or (row == 1 and col == 8) + or (row == 2 and col == 7) + or (row == 3 and col == 6) + ): mult = 0.5 else: mult = 1.0 - rch3_period_array.append(((0, row, col), 'rch_3', mult)) + rch3_period_array.append(((0, row, col), "rch_3", mult)) rch3_period[0] = rch3_period_array - rch3_package = ModflowGwfrch(model, filename='AdvGW_tidal_3.rch', - pname='rch_3', fixed_cell=True, - auxiliary='MULTIPLIER', - auxmultname='MULTIPLIER', - print_input=True, print_flows=True, - save_flows=True, - maxbound=54, - stress_period_data=rch3_period) - ts_data = [(0.0, 0.0017), (1.0, 0.0020), (11.0, 0.0017), - (21.0, 0.0018), (31.0, 0.0020)] - rch3_package.ts.initialize(timeseries=ts_data, - filename='recharge_rates_3.ts', - time_series_namerecord='rch_3', - interpolation_methodrecord='linear') + rch3_package = ModflowGwfrch( + model, + filename="AdvGW_tidal_3.rch", + pname="rch_3", + fixed_cell=True, + auxiliary="MULTIPLIER", + auxmultname="MULTIPLIER", + print_input=True, + print_flows=True, + save_flows=True, + maxbound=54, + stress_period_data=rch3_period, + ) + ts_data = [ + (0.0, 0.0017), + (1.0, 0.0020), + (11.0, 0.0017), + (21.0, 0.0018), + (31.0, 0.0020), + ] + rch3_package.ts.initialize( + timeseries=ts_data, + filename="recharge_rates_3.ts", + time_series_namerecord="rch_3", + interpolation_methodrecord="linear", + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -961,60 +1453,68 @@ def test005_advgw_tidal(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'AdvGW_tidal.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "AdvGW_tidal.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # test rename all - model.rename_all_packages('new_name') - assert model.name_file.filename == 'new_name.nam' + model.rename_all_packages("new_name") + assert model.name_file.filename == "new_name.nam" package_type_dict = {} for package in model.packagelist: if not package.package_type in package_type_dict: - assert package.filename == 'new_name.{}'.format(package.package_type) + assert package.filename == "new_name.{}".format( + package.package_type + ) package_type_dict[package.package_type] = 1 sim.write_simulation() - name_file = os.path.join(run_folder, 'new_name.nam') + name_file = os.path.join(run_folder, "new_name.nam") assert os.path.exists(name_file) - dis_file = os.path.join(run_folder, 'new_name.dis') + dis_file = os.path.join(run_folder, "new_name.dis") assert os.path.exists(dis_file) - sim.rename_all_packages('all_files_same_name') + sim.rename_all_packages("all_files_same_name") package_type_dict = {} for package in model.packagelist: if not package.package_type in package_type_dict: - assert package.filename == \ - 'all_files_same_name.{}'.format(package.package_type) + assert package.filename == "all_files_same_name.{}".format( + package.package_type + ) package_type_dict[package.package_type] = 1 - assert sim._tdis_file.filename == 'all_files_same_name.tdis' + assert sim._tdis_file.filename == "all_files_same_name.tdis" for ims_file in sim._ims_files.values(): - assert ims_file.filename == 'all_files_same_name.ims' + assert ims_file.filename == "all_files_same_name.ims" sim.write_simulation() - name_file = os.path.join(run_folder, 'all_files_same_name.nam') + name_file = os.path.join(run_folder, "all_files_same_name.nam") assert os.path.exists(name_file) - dis_file = os.path.join(run_folder, 'all_files_same_name.dis') + dis_file = os.path.join(run_folder, "all_files_same_name.dis") assert os.path.exists(dis_file) - tdis_file = os.path.join(run_folder, 'all_files_same_name.tdis') + tdis_file = os.path.join(run_folder, "all_files_same_name.tdis") assert os.path.exists(tdis_file) # load simulation - sim_load = MFSimulation.load(sim.name, 'mf6', exe_name, - sim.simulation_data.mfpath.get_sim_path(), - verbosity_level=0) + sim_load = MFSimulation.load( + sim.name, + "mf6", + exe_name, + sim.simulation_data.mfpath.get_sim_path(), + verbosity_level=0, + ) model = sim_load.get_model() # confirm ghb obs data has two blocks with correct file names - ghb = model.get_package('ghb') + ghb = model.get_package("ghb") obs = ghb.obs obs_data = obs.continuous.get_data() found_flows = False found_obs = False for key, value in obs_data.items(): - if key.lower() == 'ghb_flows.csv': + if key.lower() == "ghb_flows.csv": # there should be only one assert not found_flows found_flows = True - if key.lower() == 'ghb_obs.csv': + if key.lower() == "ghb_obs.csv": # there should be only one assert not found_obs found_obs = True @@ -1025,50 +1525,69 @@ def test005_advgw_tidal(): # check packages chk = sim.check() - summary = '.'.join(chk[0].summary_array.desc) - assert summary == '' + summary = ".".join(chk[0].summary_array.desc) + assert summary == "" return def test004_bcfss(): # init paths - test_ex_name = 'test004_bcfss' - model_name = 'bcf2ss' + test_ex_name = "test004_bcfss" + model_name = "bcf2ss" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'bcf2ss.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "bcf2ss.hds") # create simulation - sim = MFSimulation(sim_name=model_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=model_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) tdis_rc = [(1.0, 1, 1.0), (1.0, 1, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=2, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='ALL', - csv_output_filerecord='bcf2ss.ims.csv', - complexity='SIMPLE', - outer_hclose=0.000001, outer_maximum=500, - under_relaxation='NONE', inner_maximum=100, - inner_hclose=0.000001, rcloserecord=0.001, - linear_acceleration='CG', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.97) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=2, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="ALL", + csv_output_filerecord="bcf2ss.ims.csv", + complexity="SIMPLE", + outer_hclose=0.000001, + outer_maximum=500, + under_relaxation="NONE", + inner_maximum=100, + inner_hclose=0.000001, + rcloserecord=0.001, + linear_acceleration="CG", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.97, + ) sim.register_ims_package(ims_package, [model.name]) - dis_package = ModflowGwfdis(model, nlay=2, nrow=10, ncol=15, delr=500.0, - delc=500.0, - top=150.0, botm=[50.0, -50.0], - filename='{}.dis'.format(model_name)) - ic_package = ModflowGwfic(model, strt=0.0, - filename='{}.ic'.format(model_name)) + dis_package = ModflowGwfdis( + model, + nlay=2, + nrow=10, + ncol=15, + delr=500.0, + delc=500.0, + top=150.0, + botm=[50.0, -50.0], + filename="{}.dis".format(model_name), + ) + ic_package = ModflowGwfic( + model, strt=0.0, filename="{}.ic".format(model_name) + ) wetdry_data = [] for row in range(0, 10): if row == 2 or row == 7: @@ -1079,47 +1598,68 @@ def test004_bcfss(): for row in range(0, 10): for col in range(0, 15): wetdry_data.append(0.0) - npf_package = ModflowGwfnpf(model, rewet_record=[ - ('WETFCT', 1.0, 'IWETIT', 1, 'IHDWET', 0)], - save_flows=True, icelltype=[1, 0], - wetdry=wetdry_data, k=[10.0, 5.0], - k33=0.1) - oc_package = ModflowGwfoc(model, budget_filerecord='bcf2ss.cbb', - head_filerecord='bcf2ss.hds', - headprintrecord=[('COLUMNS', 15, 'WIDTH', 12, - 'DIGITS', 2, 'GENERAL')], - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) + npf_package = ModflowGwfnpf( + model, + rewet_record=[("WETFCT", 1.0, "IWETIT", 1, "IHDWET", 0)], + save_flows=True, + icelltype=[1, 0], + wetdry=wetdry_data, + k=[10.0, 5.0], + k33=0.1, + ) + oc_package = ModflowGwfoc( + model, + budget_filerecord="bcf2ss.cbb", + head_filerecord="bcf2ss.hds", + headprintrecord=[("COLUMNS", 15, "WIDTH", 12, "DIGITS", 2, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) aux = {0: [[50.0], [1.3]], 1: [[200.0], [1.5]]} # aux = {0: [[100.0], [2.3]]} - rch_package = ModflowGwfrcha(model, readasarrays=True, save_flows=True, - auxiliary=[('var1', 'var2')], - recharge={0: 0.004}, aux=aux) # *** test if aux works *** + rch_package = ModflowGwfrcha( + model, + readasarrays=True, + save_flows=True, + auxiliary=[("var1", "var2")], + recharge={0: 0.004}, + aux=aux, + ) # *** test if aux works *** # aux tests aux_out = rch_package.aux.get_data() - assert(aux_out[0][0][0,0] == 50.) - assert(aux_out[0][1][0,0] == 1.3) - assert(aux_out[1][0][0,0] == 200.0) - assert(aux_out[1][1][0,0] == 1.5) + assert aux_out[0][0][0, 0] == 50.0 + assert aux_out[0][1][0, 0] == 1.3 + assert aux_out[1][0][0, 0] == 200.0 + assert aux_out[1][1][0, 0] == 1.5 riv_period = {} riv_period_array = [] for row in range(0, 10): riv_period_array.append(((1, row, 14), 0.0, 10000.0, -5.0)) riv_period[0] = riv_period_array - riv_package = ModflowGwfriv(model, save_flows='bcf2ss.cbb', maxbound=10, - stress_period_data=riv_period) + riv_package = ModflowGwfriv( + model, + save_flows="bcf2ss.cbb", + maxbound=10, + stress_period_data=riv_period, + ) wel_period = {} - stress_period_data = [((1, 2, 3), -35000.0, 1, 2, 3), - ((1, 7, 3), -35000.0, 4, 5, 6)] + stress_period_data = [ + ((1, 2, 3), -35000.0, 1, 2, 3), + ((1, 7, 3), -35000.0, 4, 5, 6), + ] wel_period[1] = stress_period_data - wel_package = ModflowGwfwel(model, print_input=True, print_flows=True, - save_flows=True, - auxiliary=[('var1', 'var2', 'var3')], - maxbound=2, - stress_period_data=wel_period) + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + save_flows=True, + auxiliary=[("var1", "var2", "var3")], + maxbound=2, + stress_period_data=wel_period, + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -1134,10 +1674,11 @@ def test004_bcfss(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'bcf2ss.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "bcf2ss.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # clean up sim.delete_output_files() @@ -1147,76 +1688,125 @@ def test004_bcfss(): def test035_fhb(): # init paths - test_ex_name = 'test035_fhb' - model_name = 'fhb2015' + test_ex_name = "test035_fhb" + model_name = "fhb2015" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, - 'fhb2015_fhb.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join( + expected_output_folder, "fhb2015_fhb.hds" + ) # create simulation - sim = MFSimulation(sim_name=model_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=model_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) tdis_rc = [(400.0, 10, 1.0), (200.0, 4, 1.0), (400.0, 6, 1.1)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=3, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='SUMMARY', complexity='SIMPLE', - outer_hclose=0.001, - outer_maximum=120, under_relaxation='NONE', - inner_maximum=100, inner_hclose=0.0001, - rcloserecord=0.1, linear_acceleration='CG', - preconditioner_levels=7, - preconditioner_drop_tolerance=0.001, - number_orthogonalizations=2) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=3, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + complexity="SIMPLE", + outer_hclose=0.001, + outer_maximum=120, + under_relaxation="NONE", + inner_maximum=100, + inner_hclose=0.0001, + rcloserecord=0.1, + linear_acceleration="CG", + preconditioner_levels=7, + preconditioner_drop_tolerance=0.001, + number_orthogonalizations=2, + ) sim.register_ims_package(ims_package, [model.name]) - dis_package = ModflowGwfdis(model, length_units='UNDEFINED', nlay=1, - nrow=3, ncol=10, delr=1000.0, - delc=1000.0, top=50.0, botm=-200.0, - filename='{}.dis'.format(model_name)) - ic_package = ModflowGwfic(model, strt=0.0, - filename='{}.ic'.format(model_name)) - npf_package = ModflowGwfnpf(model, perched=True, icelltype=0, k=20.0, - k33=1.0) - oc_package = ModflowGwfoc(model, head_filerecord='fhb2015_fhb.hds', - headprintrecord=[('COLUMNS', 20, 'WIDTH', 5, - 'DIGITS', 2, 'FIXED')], - saverecord={0: [('HEAD', 'ALL')], - 2: [('HEAD', 'ALL')]}, - printrecord={ - 0: [('HEAD', 'ALL'), ('BUDGET', 'ALL')], - 2: [('HEAD', 'ALL'), ('BUDGET', 'ALL')]}) - sto_package = ModflowGwfsto(model, storagecoefficient=True, iconvert=0, - ss=0.01, sy=0.0) + dis_package = ModflowGwfdis( + model, + length_units="UNDEFINED", + nlay=1, + nrow=3, + ncol=10, + delr=1000.0, + delc=1000.0, + top=50.0, + botm=-200.0, + filename="{}.dis".format(model_name), + ) + ic_package = ModflowGwfic( + model, strt=0.0, filename="{}.ic".format(model_name) + ) + npf_package = ModflowGwfnpf( + model, perched=True, icelltype=0, k=20.0, k33=1.0 + ) + oc_package = ModflowGwfoc( + model, + head_filerecord="fhb2015_fhb.hds", + headprintrecord=[("COLUMNS", 20, "WIDTH", 5, "DIGITS", 2, "FIXED")], + saverecord={0: [("HEAD", "ALL")], 2: [("HEAD", "ALL")]}, + printrecord={ + 0: [("HEAD", "ALL"), ("BUDGET", "ALL")], + 2: [("HEAD", "ALL"), ("BUDGET", "ALL")], + }, + ) + sto_package = ModflowGwfsto( + model, storagecoefficient=True, iconvert=0, ss=0.01, sy=0.0 + ) time = model.modeltime - assert (time.steady_state[0] == False and time.steady_state[1] == False - and time.steady_state[2] == False) - wel_period = {0: [((0, 1, 0), 'flow')]} - wel_package = ModflowGwfwel(model, print_input=True, print_flows=True, - save_flows=True, - maxbound=1, stress_period_data=wel_period) - well_ts = [(0.0, 2000.0), (307.0, 6000.0), (791.0, 5000.0), - (1000.0, 9000.0)] - wel_package.ts.initialize(filename='fhb_flow.ts', timeseries=well_ts, - time_series_namerecord='flow', - interpolation_methodrecord='linear') + assert ( + time.steady_state[0] == False + and time.steady_state[1] == False + and time.steady_state[2] == False + ) + wel_period = {0: [((0, 1, 0), "flow")]} + wel_package = ModflowGwfwel( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=1, + stress_period_data=wel_period, + ) + well_ts = [ + (0.0, 2000.0), + (307.0, 6000.0), + (791.0, 5000.0), + (1000.0, 9000.0), + ] + wel_package.ts.initialize( + filename="fhb_flow.ts", + timeseries=well_ts, + time_series_namerecord="flow", + interpolation_methodrecord="linear", + ) chd_period = { - 0: [((0, 0, 9), 'head'), ((0, 1, 9), 'head'), ((0, 2, 9), 'head')]} - chd_package = ModflowGwfchd(model, print_input=True, print_flows=True, - save_flows=True, maxbound=3, - stress_period_data=chd_period) + 0: [((0, 0, 9), "head"), ((0, 1, 9), "head"), ((0, 2, 9), "head")] + } + chd_package = ModflowGwfchd( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=3, + stress_period_data=chd_period, + ) chd_ts = [(0.0, 0.0), (307.0, 1.0), (791.0, 5.0), (1000.0, 2.0)] - chd_package.ts.initialize(filename='fhb_head.ts', timeseries=chd_ts, - time_series_namerecord='head', - interpolation_methodrecord='linearend') + chd_package.ts.initialize( + filename="fhb_head.ts", + timeseries=chd_ts, + time_series_namerecord="head", + interpolation_methodrecord="linearend", + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -1231,10 +1821,11 @@ def test035_fhb(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'fhb2015_fhb.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "fhb2015_fhb.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # clean up sim.delete_output_files() @@ -1244,60 +1835,199 @@ def test035_fhb(): def test006_gwf3_disv(): # init paths - test_ex_name = 'test006_gwf3_disv' - model_name = 'flow' + test_ex_name = "test006_gwf3_disv" + model_name = "flow" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'flow.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "flow.hds") # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=test_ex_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) tdis_rc = [(1.0, 1, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=1, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='SUMMARY', - outer_hclose=0.00000001, - outer_maximum=1000, under_relaxation='NONE', - inner_maximum=1000, - inner_hclose=0.00000001, rcloserecord=0.01, - linear_acceleration='BICGSTAB', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.97) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=1, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + outer_hclose=0.00000001, + outer_maximum=1000, + under_relaxation="NONE", + inner_maximum=1000, + inner_hclose=0.00000001, + rcloserecord=0.01, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.97, + ) sim.register_ims_package(ims_package, [model.name]) - vertices = testutils.read_vertices(os.path.join(pth, 'vertices.txt')) - c2drecarray = testutils.read_cell2d(os.path.join(pth, 'cell2d.txt')) - disv_package = ModflowGwfdisv(model, ncpl=121, nlay=1, nvert=148, top=0.0, - botm=-100.0, idomain=1, - vertices=vertices, cell2d=c2drecarray, - filename='{}.disv'.format(model_name)) - strt_list = [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, - 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0] - ic_package = ModflowGwfic(model, strt=strt_list, - filename='{}.ic'.format(model_name)) - k = {'filename': 'k.bin', 'factor': 1.0, 'data': 1.0, 'binary': 'True'} - npf_package = ModflowGwfnpf(model, save_flows=True, icelltype=0, k=k, - k33=1.0) + vertices = testutils.read_vertices(os.path.join(pth, "vertices.txt")) + c2drecarray = testutils.read_cell2d(os.path.join(pth, "cell2d.txt")) + disv_package = ModflowGwfdisv( + model, + ncpl=121, + nlay=1, + nvert=148, + top=0.0, + botm=-100.0, + idomain=1, + vertices=vertices, + cell2d=c2drecarray, + filename="{}.disv".format(model_name), + ) + strt_list = [ + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + ] + ic_package = ModflowGwfic( + model, strt=strt_list, filename="{}.ic".format(model_name) + ) + k = {"filename": "k.bin", "factor": 1.0, "data": 1.0, "binary": "True"} + npf_package = ModflowGwfnpf( + model, save_flows=True, icelltype=0, k=k, k33=1.0 + ) k_data = npf_package.k.get_data() - assert(k_data[0,0] == 1.0) + assert k_data[0, 0] == 1.0 - oc_package = ModflowGwfoc(model, budget_filerecord='flow.cbc', - head_filerecord='flow.hds', - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) + oc_package = ModflowGwfoc( + model, + budget_filerecord="flow.cbc", + head_filerecord="flow.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) # build stress_period_data for chd package set_1 = [0, 7, 14, 18, 22, 26, 33] @@ -1307,45 +2037,58 @@ def test006_gwf3_disv(): stress_period_data.append(((0, value), 1.0)) for value in set_2: stress_period_data.append(((0, value), 0.0)) - chd_package = ModflowGwfchd(model, print_input=True, print_flows=True, - save_flows=True, maxbound=14, - stress_period_data=stress_period_data) + chd_package = ModflowGwfchd( + model, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=14, + stress_period_data=stress_period_data, + ) period_rch = {} rch_array = [] for val in range(0, 10): rch_array.append(((0, val), 0.0)) period_rch[0] = rch_array - rch_package = ModflowGwfrch(model, fixed_cell=True, maxbound=10, - stress_period_data=period_rch) - - gncrecarray = [((0, 9), (0, 40), (0, 8), 0.333333333333), - ((0, 9), (0, 42), (0, 10), 0.333333333333), - ((0, 10), (0, 43), (0, 9), 0.333333333333), - ((0, 10), (0, 45), (0, 11), 0.333333333333), - ((0, 11), (0, 46), (0, 10), 0.333333333333), - ((0, 11), (0, 48), (0, 12), 0.333333333333), - ((0, 15), (0, 40), (0, 8), 0.333333333333), - ((0, 15), (0, 58), (0, 19), 0.333333333333), - ((0, 16), (0, 48), (0, 12), 0.333333333333), - ((0, 16), (0, 66), (0, 20), 0.333333333333), - ((0, 19), (0, 67), (0, 15), 0.333333333333), - ((0, 19), (0, 85), (0, 23), 0.333333333333), - ((0, 20), (0, 75), (0, 16), 0.333333333333), - ((0, 20), (0, 93), (0, 24), 0.333333333333), - ((0, 23), (0, 94), (0, 19), 0.333333333333), - ((0, 23), (0, 112), (0, 27), 0.333333333333), - ((0, 24), (0, 102), (0, 20), 0.333333333333), - ((0, 24), (0, 120), (0, 31), 0.333333333333), - ((0, 28), (0, 112), (0, 27), 0.333333333333), - ((0, 28), (0, 114), (0, 29), 0.333333333333), - ((0, 29), (0, 115), (0, 28), 0.333333333333), - ((0, 29), (0, 117), (0, 30), 0.333333333333), - ((0, 30), (0, 118), (0, 29), 0.333333333333), - ((0, 30), (0, 120), (0, 31), 0.333333333333)] - gnc_package = ModflowGwfgnc(model, print_input=True, print_flows=True, - numgnc=24, numalphaj=1, - gncdata=gncrecarray) + rch_package = ModflowGwfrch( + model, fixed_cell=True, maxbound=10, stress_period_data=period_rch + ) + + gncrecarray = [ + ((0, 9), (0, 40), (0, 8), 0.333333333333), + ((0, 9), (0, 42), (0, 10), 0.333333333333), + ((0, 10), (0, 43), (0, 9), 0.333333333333), + ((0, 10), (0, 45), (0, 11), 0.333333333333), + ((0, 11), (0, 46), (0, 10), 0.333333333333), + ((0, 11), (0, 48), (0, 12), 0.333333333333), + ((0, 15), (0, 40), (0, 8), 0.333333333333), + ((0, 15), (0, 58), (0, 19), 0.333333333333), + ((0, 16), (0, 48), (0, 12), 0.333333333333), + ((0, 16), (0, 66), (0, 20), 0.333333333333), + ((0, 19), (0, 67), (0, 15), 0.333333333333), + ((0, 19), (0, 85), (0, 23), 0.333333333333), + ((0, 20), (0, 75), (0, 16), 0.333333333333), + ((0, 20), (0, 93), (0, 24), 0.333333333333), + ((0, 23), (0, 94), (0, 19), 0.333333333333), + ((0, 23), (0, 112), (0, 27), 0.333333333333), + ((0, 24), (0, 102), (0, 20), 0.333333333333), + ((0, 24), (0, 120), (0, 31), 0.333333333333), + ((0, 28), (0, 112), (0, 27), 0.333333333333), + ((0, 28), (0, 114), (0, 29), 0.333333333333), + ((0, 29), (0, 115), (0, 28), 0.333333333333), + ((0, 29), (0, 117), (0, 30), 0.333333333333), + ((0, 30), (0, 118), (0, 29), 0.333333333333), + ((0, 30), (0, 120), (0, 31), 0.333333333333), + ] + gnc_package = ModflowGwfgnc( + model, + print_input=True, + print_flows=True, + numgnc=24, + numalphaj=1, + gncdata=gncrecarray, + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -1359,13 +2102,14 @@ def test006_gwf3_disv(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'flow.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "flow.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # export to netcdf - temporarily disabled - #model.export(os.path.join(run_folder, "test006_gwf3.nc")) + # model.export(os.path.join(run_folder, "test006_gwf3.nc")) # export to shape file model.export(os.path.join(run_folder, "test006_gwf3.shp")) @@ -1377,82 +2121,208 @@ def test006_gwf3_disv(): def test006_2models_gnc(): # init paths - test_ex_name = 'test006_2models_gnc' - model_name_1 = 'model1' - model_name_2 = 'model2' + test_ex_name = "test006_2models_gnc" + model_name_1 = "model1" + model_name_2 = "model2" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file_1 = os.path.join(expected_output_folder, 'model1.hds') - expected_head_file_2 = os.path.join(expected_output_folder, 'model2.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file_1 = os.path.join(expected_output_folder, "model1.hds") + expected_head_file_2 = os.path.join(expected_output_folder, "model2.hds") # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=test_ex_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) tdis_rc = [(1.0, 1, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=1, - perioddata=tdis_rc) - model_1 = ModflowGwf(sim, modelname=model_name_1, - model_nam_file='{}.nam'.format(model_name_1)) - model_2 = ModflowGwf(sim, modelname=model_name_2, - model_nam_file='{}.nam'.format(model_name_2)) - ims_package = ModflowIms(sim, print_option='SUMMARY', - outer_hclose=0.00000001, - outer_maximum=1000, under_relaxation='NONE', - inner_maximum=1000, - inner_hclose=0.00000001, rcloserecord=0.01, - linear_acceleration='BICGSTAB', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.97) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=1, perioddata=tdis_rc + ) + model_1 = ModflowGwf( + sim, + modelname=model_name_1, + model_nam_file="{}.nam".format(model_name_1), + ) + model_2 = ModflowGwf( + sim, + modelname=model_name_2, + model_nam_file="{}.nam".format(model_name_2), + ) + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + outer_hclose=0.00000001, + outer_maximum=1000, + under_relaxation="NONE", + inner_maximum=1000, + inner_hclose=0.00000001, + rcloserecord=0.01, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.97, + ) sim.register_ims_package(ims_package, [model_1.name, model_2.name]) - idom = [1, 1, 1, 1, 1, 1, 1, - 1, 1, 1, 1, 1, 1, 1, - 1, 1, 0, 0, 0, 1, 1, - 1, 1, 0, 0, 0, 1, 1, - 1, 1, 0, 0, 0, 1, 1, - 1, 1, 1, 1, 1, 1, 1, - 1, 1, 1, 1, 1, 1, 1, ] - dis_package_1 = ModflowGwfdis(model_1, length_units='METERS', nlay=1, - nrow=7, ncol=7, idomain=idom, - delr=100.0, delc=100.0, top=0.0, botm=-100.0, - filename='{}.dis'.format(model_name_1)) - dis_package_2 = ModflowGwfdis(model_2, length_units='METERS', nlay=1, - nrow=9, ncol=9, delr=33.33, - delc=33.33, top=0.0, botm=-100.0, - filename='{}.dis'.format(model_name_2)) - - strt_list = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, ] - ic_package_1 = ModflowGwfic(model_1, strt=strt_list, - filename='{}.ic'.format(model_name_1)) - ic_package_2 = ModflowGwfic(model_2, strt=1.0, - filename='{}.ic'.format(model_name_2)) - npf_package_1 = ModflowGwfnpf(model_1, save_flows=True, perched=True, - icelltype=0, k=1.0, k33=1.0) - npf_package_2 = ModflowGwfnpf(model_2, save_flows=True, perched=True, - icelltype=0, k=1.0, k33=1.0) - oc_package_1 = ModflowGwfoc(model_1, budget_filerecord='model1.cbc', - head_filerecord='model1.hds', - saverecord=[('HEAD', 'ALL'), - ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), - ('BUDGET', 'ALL')]) - oc_package_2 = ModflowGwfoc(model_2, budget_filerecord='model2.cbc', - head_filerecord='model2.hds', - saverecord=[('HEAD', 'ALL'), - ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), - ('BUDGET', 'ALL')]) + idom = [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 0, + 0, + 0, + 1, + 1, + 1, + 1, + 0, + 0, + 0, + 1, + 1, + 1, + 1, + 0, + 0, + 0, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + ] + dis_package_1 = ModflowGwfdis( + model_1, + length_units="METERS", + nlay=1, + nrow=7, + ncol=7, + idomain=idom, + delr=100.0, + delc=100.0, + top=0.0, + botm=-100.0, + filename="{}.dis".format(model_name_1), + ) + dis_package_2 = ModflowGwfdis( + model_2, + length_units="METERS", + nlay=1, + nrow=9, + ncol=9, + delr=33.33, + delc=33.33, + top=0.0, + botm=-100.0, + filename="{}.dis".format(model_name_2), + ) + + strt_list = [ + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.0, + ] + ic_package_1 = ModflowGwfic( + model_1, strt=strt_list, filename="{}.ic".format(model_name_1) + ) + ic_package_2 = ModflowGwfic( + model_2, strt=1.0, filename="{}.ic".format(model_name_2) + ) + npf_package_1 = ModflowGwfnpf( + model_1, save_flows=True, perched=True, icelltype=0, k=1.0, k33=1.0 + ) + npf_package_2 = ModflowGwfnpf( + model_2, save_flows=True, perched=True, icelltype=0, k=1.0, k33=1.0 + ) + oc_package_1 = ModflowGwfoc( + model_1, + budget_filerecord="model1.cbc", + head_filerecord="model1.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + oc_package_2 = ModflowGwfoc( + model_2, + budget_filerecord="model2.cbc", + head_filerecord="model2.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) # build periodrecarray for chd package set_1 = [0, 7, 14, 18, 22, 26, 33] @@ -1462,46 +2332,78 @@ def test006_2models_gnc(): stress_period_data.append(((0, value, 0), 1.0)) for value in range(0, 7): stress_period_data.append(((0, value, 6), 0.0)) - chd_package = ModflowGwfchd(model_1, print_input=True, print_flows=True, - save_flows=True, maxbound=30, - stress_period_data=stress_period_data) - - gncrecarray = testutils.read_gncrecarray(os.path.join(pth, 'gnc.txt')) + chd_package = ModflowGwfchd( + model_1, + print_input=True, + print_flows=True, + save_flows=True, + maxbound=30, + stress_period_data=stress_period_data, + ) + + gncrecarray = testutils.read_gncrecarray(os.path.join(pth, "gnc.txt")) # test gnc delete new_gncrecarray = gncrecarray[10:] - gnc_package = ModflowGwfgnc(sim, print_input=True, print_flows=True, - numgnc=26, numalphaj=1, - gncdata=new_gncrecarray) + gnc_package = ModflowGwfgnc( + sim, + print_input=True, + print_flows=True, + numgnc=26, + numalphaj=1, + gncdata=new_gncrecarray, + ) sim.remove_package(gnc_package.package_type) - gnc_package = ModflowGwfgnc(sim, print_input=True, print_flows=True, - numgnc=36, numalphaj=1, - gncdata=gncrecarray) + gnc_package = ModflowGwfgnc( + sim, + print_input=True, + print_flows=True, + numgnc=36, + numalphaj=1, + gncdata=gncrecarray, + ) - exgrecarray = testutils.read_exchangedata(os.path.join(pth, 'exg.txt')) + exgrecarray = testutils.read_exchangedata(os.path.join(pth, "exg.txt")) # build obs dictionary - gwf_obs = {('gwfgwf_obs.csv'): [('gwf-1-3-2_1-1-1', 'flow-ja-face', - (0, 2, 1), (0, 0, 0)), - ('gwf-1-3-2_1-2-1', 'flow-ja-face', - (0, 2, 1), (0, 1, 0))]} + gwf_obs = { + ("gwfgwf_obs.csv"): [ + ("gwf-1-3-2_1-1-1", "flow-ja-face", (0, 2, 1), (0, 0, 0)), + ("gwf-1-3-2_1-2-1", "flow-ja-face", (0, 2, 1), (0, 1, 0)), + ] + } # test exg delete newexgrecarray = exgrecarray[10:] - exg_package = ModflowGwfgwf(sim, print_input=True, print_flows=True, - save_flows=True, auxiliary='testaux', - gnc_filerecord='test006_2models_gnc.gnc', - nexg=26, exchangedata=newexgrecarray, - exgtype='gwf6-gwf6', exgmnamea=model_name_1, - exgmnameb=model_name_2) + exg_package = ModflowGwfgwf( + sim, + print_input=True, + print_flows=True, + save_flows=True, + auxiliary="testaux", + gnc_filerecord="test006_2models_gnc.gnc", + nexg=26, + exchangedata=newexgrecarray, + exgtype="gwf6-gwf6", + exgmnamea=model_name_1, + exgmnameb=model_name_2, + ) sim.remove_package(exg_package.package_type) - exg_package = ModflowGwfgwf(sim, print_input=True, print_flows=True, - save_flows=True, auxiliary='testaux', - gnc_filerecord='test006_2models_gnc.gnc', - nexg=36, exchangedata=exgrecarray, - exgtype='gwf6-gwf6', exgmnamea=model_name_1, - exgmnameb=model_name_2, observations=gwf_obs) + exg_package = ModflowGwfgwf( + sim, + print_input=True, + print_flows=True, + save_flows=True, + auxiliary="testaux", + gnc_filerecord="test006_2models_gnc.gnc", + nexg=36, + exchangedata=exgrecarray, + exgtype="gwf6-gwf6", + exgmnamea=model_name_1, + exgmnameb=model_name_2, + observations=gwf_obs, + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -1516,17 +2418,19 @@ def test006_2models_gnc(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file_1) - head_new = os.path.join(run_folder, 'model1.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "model1.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file_2) - head_new = os.path.join(run_folder, 'model2.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "model2.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # clean up sim.delete_output_files() @@ -1536,58 +2440,82 @@ def test006_2models_gnc(): def test050_circle_island(): # init paths - test_ex_name = 'test050_circle_island' - model_name = 'ci' + test_ex_name = "test050_circle_island" + model_name = "ci" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'ci.output.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "ci.output.hds") # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=test_ex_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) tdis_rc = [(1.0, 1, 1.0)] - tdis_package = ModflowTdis(sim, time_units='DAYS', nper=1, - perioddata=tdis_rc) - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) - ims_package = ModflowIms(sim, print_option='SUMMARY', - outer_hclose=0.000001, - outer_maximum=500, under_relaxation='NONE', - inner_maximum=1000, - inner_hclose=0.000001, rcloserecord=0.000001, - linear_acceleration='BICGSTAB', - relaxation_factor=0.0) + tdis_package = ModflowTdis( + sim, time_units="DAYS", nper=1, perioddata=tdis_rc + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + outer_hclose=0.000001, + outer_maximum=500, + under_relaxation="NONE", + inner_maximum=1000, + inner_hclose=0.000001, + rcloserecord=0.000001, + linear_acceleration="BICGSTAB", + relaxation_factor=0.0, + ) sim.register_ims_package(ims_package, [model.name]) - vertices = testutils.read_vertices(os.path.join(pth, 'vertices.txt')) - c2drecarray = testutils.read_cell2d(os.path.join(pth, 'cell2d.txt')) - disv_package = ModflowGwfdisv(model, ncpl=5240, nlay=2, nvert=2778, - top=0.0, botm=[-20.0, -40.0], - idomain=1, vertices=vertices, - cell2d=c2drecarray, - filename='{}.disv'.format(model_name)) - ic_package = ModflowGwfic(model, strt=0.0, - filename='{}.ic'.format(model_name)) - npf_package = ModflowGwfnpf(model, save_flows=True, icelltype=0, k=10.0, - k33=0.2) - oc_package = ModflowGwfoc(model, budget_filerecord='ci.output.cbc', - head_filerecord='ci.output.hds', - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) + vertices = testutils.read_vertices(os.path.join(pth, "vertices.txt")) + c2drecarray = testutils.read_cell2d(os.path.join(pth, "cell2d.txt")) + disv_package = ModflowGwfdisv( + model, + ncpl=5240, + nlay=2, + nvert=2778, + top=0.0, + botm=[-20.0, -40.0], + idomain=1, + vertices=vertices, + cell2d=c2drecarray, + filename="{}.disv".format(model_name), + ) + ic_package = ModflowGwfic( + model, strt=0.0, filename="{}.ic".format(model_name) + ) + npf_package = ModflowGwfnpf( + model, save_flows=True, icelltype=0, k=10.0, k33=0.2 + ) + oc_package = ModflowGwfoc( + model, + budget_filerecord="ci.output.cbc", + head_filerecord="ci.output.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) stress_period_data = testutils.read_ghbrecarray( - os.path.join(pth, 'ghb.txt'), 2) - ghb_package = ModflowGwfghb(model, maxbound=3173, - stress_period_data=stress_period_data) + os.path.join(pth, "ghb.txt"), 2 + ) + ghb_package = ModflowGwfghb( + model, maxbound=3173, stress_period_data=stress_period_data + ) - rch_data = ['OPEN/CLOSE', 'rech.dat', 'FACTOR', 1.0, 'IPRN', 0] - rch_package = ModflowGwfrcha(model, readasarrays=True, - save_flows=True, recharge=rch_data) + rch_data = ["OPEN/CLOSE", "rech.dat", "FACTOR", 1.0, "IPRN", 0] + rch_package = ModflowGwfrcha( + model, readasarrays=True, save_flows=True, recharge=rch_data + ) # change folder to save simulation sim.simulation_data.mfpath.set_sim_path(run_folder) @@ -1602,10 +2530,11 @@ def test050_circle_island(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'ci.output.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "ci.output.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # clean up sim.delete_output_files() @@ -1615,167 +2544,239 @@ def test050_circle_island(): def test028_sfr(): # init paths - test_ex_name = 'test028_sfr' - model_name = 'test1tr' + test_ex_name = "test028_sfr" + model_name = "test1tr" - pth = os.path.join('..', 'examples', 'data', 'mf6', 'create_tests', - test_ex_name) + pth = os.path.join( + "..", "examples", "data", "mf6", "create_tests", test_ex_name + ) run_folder = os.path.join(cpth, test_ex_name) if not os.path.isdir(run_folder): os.makedirs(run_folder) - expected_output_folder = os.path.join(pth, 'expected_output') - expected_head_file = os.path.join(expected_output_folder, 'test1tr.hds') + expected_output_folder = os.path.join(pth, "expected_output") + expected_head_file = os.path.join(expected_output_folder, "test1tr.hds") # create simulation - sim = MFSimulation(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=pth) + sim = MFSimulation( + sim_name=test_ex_name, version="mf6", exe_name=exe_name, sim_ws=pth + ) sim.name_file.continue_.set_data(True) tdis_rc = [(1577889000, 50, 1.1), (1577889000, 50, 1.1)] - tdis_package = ModflowTdis(sim, time_units='SECONDS', nper=2, - perioddata=tdis_rc, filename='simulation.tdis') - model = ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) + tdis_package = ModflowTdis( + sim, + time_units="SECONDS", + nper=2, + perioddata=tdis_rc, + filename="simulation.tdis", + ) + model = ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) model.name_file.save_flows.set_data(True) - ims_package = ModflowIms(sim, print_option='SUMMARY', outer_hclose=0.00001, - outer_maximum=100, under_relaxation='DBD', - under_relaxation_theta=0.85, - under_relaxation_kappa=0.0001, - under_relaxation_gamma=0.0, - under_relaxation_momentum=0.1, - backtracking_number=0, backtracking_tolerance=1.1, - backtracking_reduction_factor=0.7, - backtracking_residual_limit=1.0, - inner_hclose=0.00001, rcloserecord=0.1, - inner_maximum=100, linear_acceleration='CG', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.99, - filename='model.ims') + ims_package = ModflowIms( + sim, + print_option="SUMMARY", + outer_hclose=0.00001, + outer_maximum=100, + under_relaxation="DBD", + under_relaxation_theta=0.85, + under_relaxation_kappa=0.0001, + under_relaxation_gamma=0.0, + under_relaxation_momentum=0.1, + backtracking_number=0, + backtracking_tolerance=1.1, + backtracking_reduction_factor=0.7, + backtracking_residual_limit=1.0, + inner_hclose=0.00001, + rcloserecord=0.1, + inner_maximum=100, + linear_acceleration="CG", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.99, + filename="model.ims", + ) sim.register_ims_package(ims_package, [model.name]) - top = testutils.read_std_array(os.path.join(pth, 'top.txt'), 'float') - botm = testutils.read_std_array(os.path.join(pth, 'botm.txt'), 'float') - idomain = testutils.read_std_array(os.path.join(pth, 'idomain.txt'), 'int') - dis_package = ModflowGwfdis(model, length_units='FEET', nlay=1, nrow=15, - ncol=10, delr=5000.0, delc=5000.0, - top=top, botm=botm, idomain=idomain, - filename='{}.dis'.format(model_name)) - strt = testutils.read_std_array(os.path.join(pth, 'strt.txt'), 'float') - strt_int = ['internal', 'factor', 1.0, 'iprn', 0, strt] - ic_package = ModflowGwfic(model, strt=strt_int, - filename='{}.ic'.format(model_name)) - - k_vals = testutils.read_std_array(os.path.join(pth, 'k.txt'), 'float') - k = ['internal', 'factor', 3.000E-03, 'iprn', 0, k_vals] + top = testutils.read_std_array(os.path.join(pth, "top.txt"), "float") + botm = testutils.read_std_array(os.path.join(pth, "botm.txt"), "float") + idomain = testutils.read_std_array(os.path.join(pth, "idomain.txt"), "int") + dis_package = ModflowGwfdis( + model, + length_units="FEET", + nlay=1, + nrow=15, + ncol=10, + delr=5000.0, + delc=5000.0, + top=top, + botm=botm, + idomain=idomain, + filename="{}.dis".format(model_name), + ) + strt = testutils.read_std_array(os.path.join(pth, "strt.txt"), "float") + strt_int = ["internal", "factor", 1.0, "iprn", 0, strt] + ic_package = ModflowGwfic( + model, strt=strt_int, filename="{}.ic".format(model_name) + ) + + k_vals = testutils.read_std_array(os.path.join(pth, "k.txt"), "float") + k = ["internal", "factor", 3.000e-03, "iprn", 0, k_vals] npf_package = ModflowGwfnpf(model, icelltype=1, k=k, k33=1.0) - npf_package.k.factor = 2.000E-04 + npf_package.k.factor = 2.000e-04 - oc_package = ModflowGwfoc(model, budget_filerecord='test1tr.cbc', - head_filerecord='test1tr.hds', - saverecord={0: [('HEAD', 'FREQUENCY', 5), - ('BUDGET', 'FREQUENCY', 5)]}, - printrecord={0: [('HEAD', 'FREQUENCY', 5), - ('BUDGET', 'FREQUENCY', 5)]}) + oc_package = ModflowGwfoc( + model, + budget_filerecord="test1tr.cbc", + head_filerecord="test1tr.hds", + saverecord={0: [("HEAD", "FREQUENCY", 5), ("BUDGET", "FREQUENCY", 5)]}, + printrecord={ + 0: [("HEAD", "FREQUENCY", 5), ("BUDGET", "FREQUENCY", 5)] + }, + ) - sy_vals = testutils.read_std_array(os.path.join(pth, 'sy.txt'), 'float') - sy = {'factor': 0.2, 'iprn': 0, 'data': sy_vals} - sto_package = ModflowGwfsto(model, iconvert=1, ss=1.0E-6, sy=sy) + sy_vals = testutils.read_std_array(os.path.join(pth, "sy.txt"), "float") + sy = {"factor": 0.2, "iprn": 0, "data": sy_vals} + sto_package = ModflowGwfsto(model, iconvert=1, ss=1.0e-6, sy=sy) - surf = testutils.read_std_array(os.path.join(pth, 'surface.txt'), 'float') - surf_data = ['internal', 'factor', 1.0, 'iprn', -1, surf] + surf = testutils.read_std_array(os.path.join(pth, "surface.txt"), "float") + surf_data = ["internal", "factor", 1.0, "iprn", -1, surf] # build time array series - tas = {0.0: 9.5E-08, 6.0E09: 9.5E-08, - 'filename': 'test028_sfr.evt.tas', - 'time_series_namerecord': 'evtarray_1', - 'interpolation_methodrecord': 'LINEAR'} - - evt_package = ModflowGwfevta(model, readasarrays=True, timearrayseries=tas, - surface=surf_data, depth=15.0, - rate='TIMEARRAYSERIES evtarray_1', - filename='test1tr.evt') + tas = { + 0.0: 9.5e-08, + 6.0e09: 9.5e-08, + "filename": "test028_sfr.evt.tas", + "time_series_namerecord": "evtarray_1", + "interpolation_methodrecord": "LINEAR", + } + + evt_package = ModflowGwfevta( + model, + readasarrays=True, + timearrayseries=tas, + surface=surf_data, + depth=15.0, + rate="TIMEARRAYSERIES evtarray_1", + filename="test1tr.evt", + ) # attach obs package to evt - obs_dict = {'test028_sfr.evt.csv': [('obs-1', 'EVT', (0, 1, 5)), - ('obs-2', 'EVT', (0, 2, 3))]} - evt_package.obs.initialize(filename='test028_sfr.evt.obs', print_input=True, - continuous=obs_dict) + obs_dict = { + "test028_sfr.evt.csv": [ + ("obs-1", "EVT", (0, 1, 5)), + ("obs-2", "EVT", (0, 2, 3)), + ] + } + evt_package.obs.initialize( + filename="test028_sfr.evt.obs", print_input=True, continuous=obs_dict + ) stress_period_data = { - 0: [((0, 12, 0), 988.0, 0.038), ((0, 13, 8), 1045.0, 0.038)]} - ghb_package = ModflowGwfghb(model, maxbound=2, - stress_period_data=stress_period_data) + 0: [((0, 12, 0), 988.0, 0.038), ((0, 13, 8), 1045.0, 0.038)] + } + ghb_package = ModflowGwfghb( + model, maxbound=2, stress_period_data=stress_period_data + ) - rch = testutils.read_std_array(os.path.join(pth, 'recharge.txt'), 'float') + rch = testutils.read_std_array(os.path.join(pth, "recharge.txt"), "float") # test empty rch_data = ModflowGwfrcha.recharge.empty(model) - rch_data[0]['data'] = rch - rch_data[0]['factor'] = 5.000E-10 - rch_data[0]['iprn'] = -1 - rch_package = ModflowGwfrcha(model, readasarrays=True, recharge=rch_data, - filename='test1tr.rch') - - sfr_rec = testutils.read_sfr_rec(os.path.join(pth, 'sfr_rec.txt'), 3) + rch_data[0]["data"] = rch + rch_data[0]["factor"] = 5.000e-10 + rch_data[0]["iprn"] = -1 + rch_package = ModflowGwfrcha( + model, readasarrays=True, recharge=rch_data, filename="test1tr.rch" + ) + + sfr_rec = testutils.read_sfr_rec(os.path.join(pth, "sfr_rec.txt"), 3) reach_con_rec = testutils.read_reach_con_rec( - os.path.join(pth, 'sfr_reach_con_rec.txt')) + os.path.join(pth, "sfr_reach_con_rec.txt") + ) reach_div_rec = testutils.read_reach_div_rec( - os.path.join(pth, 'sfr_reach_div_rec.txt')) + os.path.join(pth, "sfr_reach_div_rec.txt") + ) reach_per_rec = testutils.read_reach_per_rec( - os.path.join(pth, 'sfr_reach_per_rec.txt')) + os.path.join(pth, "sfr_reach_per_rec.txt") + ) # test zero based indexes reach_con_rec[0] = (0, -0.0) - sfr_package = ModflowGwfsfr(model, unit_conversion=1.486, - stage_filerecord='test1tr.sfr.stage.bin', - budget_filerecord='test1tr.sfr.cbc', - nreaches=36, packagedata=sfr_rec, - connectiondata=reach_con_rec, - diversions=reach_div_rec, - perioddata={0: reach_per_rec}) - assert (sfr_package.connectiondata.get_data()[0][1] == -0.0) - assert (sfr_package.connectiondata.get_data()[1][1] == 0.0) - assert (sfr_package.connectiondata.get_data()[2][1] == 1.0) - assert (sfr_package.packagedata.get_data()[1][1].lower() == 'none') + sfr_package = ModflowGwfsfr( + model, + unit_conversion=1.486, + stage_filerecord="test1tr.sfr.stage.bin", + budget_filerecord="test1tr.sfr.cbc", + nreaches=36, + packagedata=sfr_rec, + connectiondata=reach_con_rec, + diversions=reach_div_rec, + perioddata={0: reach_per_rec}, + ) + assert sfr_package.connectiondata.get_data()[0][1] == -0.0 + assert sfr_package.connectiondata.get_data()[1][1] == 0.0 + assert sfr_package.connectiondata.get_data()[2][1] == 1.0 + assert sfr_package.packagedata.get_data()[1][1].lower() == "none" sim.simulation_data.mfpath.set_sim_path(run_folder) sim.write_simulation() - sim.load(sim_name=test_ex_name, version='mf6', exe_name=exe_name, - sim_ws=run_folder) + sim.load( + sim_name=test_ex_name, + version="mf6", + exe_name=exe_name, + sim_ws=run_folder, + ) model = sim.get_model(model_name) - sfr_package = model.get_package('sfr') + sfr_package = model.get_package("sfr") # sfr_package.set_all_data_external() - assert (sfr_package.connectiondata.get_data()[0][1] == -0.0) - assert (sfr_package.connectiondata.get_data()[1][1] == 0.0) - assert (sfr_package.connectiondata.get_data()[2][1] == 1.0) + assert sfr_package.connectiondata.get_data()[0][1] == -0.0 + assert sfr_package.connectiondata.get_data()[1][1] == 0.0 + assert sfr_package.connectiondata.get_data()[2][1] == 1.0 pdata = sfr_package.packagedata.get_data() - assert (sfr_package.packagedata.get_data()[1][1].lower() == 'none') + assert sfr_package.packagedata.get_data()[1][1].lower() == "none" # undo zero based test and move on model.remove_package(sfr_package.package_type) reach_con_rec = testutils.read_reach_con_rec( - os.path.join(pth, 'sfr_reach_con_rec.txt')) + os.path.join(pth, "sfr_reach_con_rec.txt") + ) # set sfr settings back to expected package data rec_line = (sfr_rec[1][0], (0, 1, 1)) + sfr_rec[1][2:] sfr_rec[1] = rec_line - sfr_package = ModflowGwfsfr(model, unit_conversion=1.486, - stage_filerecord='test1tr.sfr.stage.bin', - budget_filerecord='test1tr.sfr.cbc', - nreaches=36, packagedata=sfr_rec, - connectiondata=reach_con_rec, - diversions=reach_div_rec, - perioddata={0: reach_per_rec}) - - obs_data_1 = testutils.read_obs(os.path.join(pth, 'sfr_obs_1.txt')) - obs_data_2 = testutils.read_obs(os.path.join(pth, 'sfr_obs_2.txt')) - obs_data_3 = testutils.read_obs(os.path.join(pth, 'sfr_obs_3.txt')) - obs_data = {'test1tr.sfr.csv': obs_data_1, - 'test1tr.sfr.qaq.csv': obs_data_2, - 'test1tr.sfr.flow.csv': obs_data_3} - sfr_package.obs.initialize(filename='test1tr.sfr.obs', digits=10, - print_input=True, continuous=obs_data) - - wells = testutils.read_wells(os.path.join(pth, 'well.txt')) - wel_package = ModflowGwfwel(model, boundnames=True, maxbound=10, - stress_period_data={0: wells, 1: [()]}) + sfr_package = ModflowGwfsfr( + model, + unit_conversion=1.486, + stage_filerecord="test1tr.sfr.stage.bin", + budget_filerecord="test1tr.sfr.cbc", + nreaches=36, + packagedata=sfr_rec, + connectiondata=reach_con_rec, + diversions=reach_div_rec, + perioddata={0: reach_per_rec}, + ) + + obs_data_1 = testutils.read_obs(os.path.join(pth, "sfr_obs_1.txt")) + obs_data_2 = testutils.read_obs(os.path.join(pth, "sfr_obs_2.txt")) + obs_data_3 = testutils.read_obs(os.path.join(pth, "sfr_obs_3.txt")) + obs_data = { + "test1tr.sfr.csv": obs_data_1, + "test1tr.sfr.qaq.csv": obs_data_2, + "test1tr.sfr.flow.csv": obs_data_3, + } + sfr_package.obs.initialize( + filename="test1tr.sfr.obs", + digits=10, + print_input=True, + continuous=obs_data, + ) + + wells = testutils.read_wells(os.path.join(pth, "well.txt")) + wel_package = ModflowGwfwel( + model, + boundnames=True, + maxbound=10, + stress_period_data={0: wells, 1: [()]}, + ) # write simulation to new location sim.write_simulation() @@ -1786,10 +2787,11 @@ def test028_sfr(): # compare output to expected results head_file = os.path.join(os.getcwd(), expected_head_file) - head_new = os.path.join(run_folder, 'test1tr.hds') - outfile = os.path.join(run_folder, 'head_compare.dat') - assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, - outfile=outfile) + head_new = os.path.join(run_folder, "test1tr.hds") + outfile = os.path.join(run_folder, "head_compare.dat") + assert pymake.compare_heads( + None, None, files1=head_file, files2=head_new, outfile=outfile + ) # clean up sim.delete_output_files() @@ -1797,7 +2799,7 @@ def test028_sfr(): return -if __name__ == '__main__': +if __name__ == "__main__": np001() np002() test004_bcfss() From 5927f544984bc4a315f08dd74bd7ca8f9a1b7415 Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Thu, 27 Aug 2020 08:08:48 -0700 Subject: [PATCH 3/6] feat(cellid not as tuple): tests --- autotest/t505_test.py | 58 +++++++++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 24 deletions(-) diff --git a/autotest/t505_test.py b/autotest/t505_test.py index a90dc882a5..e81c5793dc 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -2056,30 +2056,30 @@ def test006_gwf3_disv(): ) gncrecarray = [ - ((0, 9), (0, 40), (0, 8), 0.333333333333), - ((0, 9), (0, 42), (0, 10), 0.333333333333), - ((0, 10), (0, 43), (0, 9), 0.333333333333), - ((0, 10), (0, 45), (0, 11), 0.333333333333), - ((0, 11), (0, 46), (0, 10), 0.333333333333), - ((0, 11), (0, 48), (0, 12), 0.333333333333), - ((0, 15), (0, 40), (0, 8), 0.333333333333), - ((0, 15), (0, 58), (0, 19), 0.333333333333), - ((0, 16), (0, 48), (0, 12), 0.333333333333), - ((0, 16), (0, 66), (0, 20), 0.333333333333), - ((0, 19), (0, 67), (0, 15), 0.333333333333), - ((0, 19), (0, 85), (0, 23), 0.333333333333), - ((0, 20), (0, 75), (0, 16), 0.333333333333), - ((0, 20), (0, 93), (0, 24), 0.333333333333), - ((0, 23), (0, 94), (0, 19), 0.333333333333), - ((0, 23), (0, 112), (0, 27), 0.333333333333), - ((0, 24), (0, 102), (0, 20), 0.333333333333), - ((0, 24), (0, 120), (0, 31), 0.333333333333), - ((0, 28), (0, 112), (0, 27), 0.333333333333), - ((0, 28), (0, 114), (0, 29), 0.333333333333), - ((0, 29), (0, 115), (0, 28), 0.333333333333), - ((0, 29), (0, 117), (0, 30), 0.333333333333), - ((0, 30), (0, 118), (0, 29), 0.333333333333), - ((0, 30), (0, 120), (0, 31), 0.333333333333), + (0, 9, 0, 40, 0, 8, 0.333333333333), + (0, 9, 0, 42, 0, 10, 0.333333333333), + (0, 10, 0, 43, 0, 9, 0.333333333333), + (0, 10, 0, 45, 0, 11, 0.333333333333), + (0, 11, 0, 46, 0, 10, 0.333333333333), + (0, 11, 0, 48, 0, 12, 0.333333333333), + (0, 15, 0, 40, 0, 8, 0.333333333333), + (0, 15, 0, 58, 0, 19, 0.333333333333), + (0, 16, 0, 48, 0, 12, 0.333333333333), + (0, 16, 0, 66, 0, 20, 0.333333333333), + (0, 19, 0, 67, 0, 15, 0.333333333333), + (0, 19, 0, 85, 0, 23, 0.333333333333), + (0, 20, 0, 75, 0, 16, 0.333333333333), + (0, 20, 0, 93, 0, 24, 0.333333333333), + (0, 23, 0, 94, 0, 19, 0.333333333333), + (0, 23, 0, 112, 0, 27, 0.333333333333), + (0, 24, 0, 102, 0, 20, 0.333333333333), + (0, 24, 0, 120, 0, 31, 0.333333333333), + (0, 28, 0, 112, 0, 27, 0.333333333333), + (0, 28, 0, 114, 0, 29, 0.333333333333), + (0, 29, 0, 115, 0, 28, 0.333333333333), + (0, 29, 0, 117, 0, 30, 0.333333333333), + (0, 30, 0, 118, 0, 29, 0.333333333333), + (0, 30, 0, 120, 0, 31, 0.333333333333), ] gnc_package = ModflowGwfgnc( model, @@ -2698,6 +2698,16 @@ def test028_sfr(): reach_per_rec = testutils.read_reach_per_rec( os.path.join(pth, "sfr_reach_per_rec.txt") ) + + # test trying to get empty perioddata + # should fail because data is jagged + try: + pd = ModflowGwfsfr.perioddata.empty(model) + success = True + except MFDataException: + success = False + assert not success + # test zero based indexes reach_con_rec[0] = (0, -0.0) sfr_package = ModflowGwfsfr( From e7ff1c8d07e82c098a3fd922a125965c176de016 Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Thu, 27 Aug 2020 09:31:28 -0700 Subject: [PATCH 4/6] fix(line spacing): black not functioning the same on my desktop as the server. trying "black --line-spacing 78" since 79 appears to not be cutting off all the lines at 79. --- flopy/mf6/data/mfdatautil.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 6aab1db991..0077996871 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -169,7 +169,10 @@ def to_string( 'Cellid "{}" contains {} integer(s). Expected a' " cellid containing {} integer(s) for grid type" " {}.".format( - val, len(val), cellid_size, str(model_grid.grid_type()) + val, + len(val), + cellid_size, + str(model_grid.grid_type()), ) ) type_, value_, traceback_ = sys.exc_info() @@ -731,7 +734,8 @@ def dtype( ) if aux_vars is not None: if len(aux_vars) > 0 and ( - isinstance(aux_vars[0], list) or isinstance(aux_vars[0], tuple) + isinstance(aux_vars[0], list) + or isinstance(aux_vars[0], tuple) ): aux_vars = aux_vars[0] for aux_var in aux_vars: From 4826cdd621a0c274eb4595a73e8841017f53f72d Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Thu, 27 Aug 2020 10:44:02 -0700 Subject: [PATCH 5/6] fix(style issues) --- flopy/mf6/data/mfdatastorage.py | 105 ++++++++++++++++---------------- flopy/mf6/data/mfdatautil.py | 12 ++-- 2 files changed, 57 insertions(+), 60 deletions(-) diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py index cf67ade0bb..5705623633 100644 --- a/flopy/mf6/data/mfdatastorage.py +++ b/flopy/mf6/data/mfdatastorage.py @@ -227,8 +227,8 @@ class DataStorage(object): data is a list containing all data for a single record in the recarray. . data structure type must be recarray append_data(data) - appends data "data" to the end of a recarray. data structure type must - be recarray + appends data "data" to the end of a recarray. data structure type + must be recarray set_data(data, layer=None, multiplier=[1.0] sets the data being stored to "data" for layer "layer", replacing all data for that layer. a multiplier can be specified. @@ -488,7 +488,8 @@ def _get_layer_header_str(self, layer): self.layer_storage[layer].data_storage_type == DataStorageType.internal_constant ): - header_list.append("constant {}".format(self.layer_storage[layer])) + lr = self.layer_storage[layer] + header_list.append("constant {}".format(lr)) else: header_list.append("internal") if ( @@ -715,8 +716,9 @@ def _access_data(self, layer, return_data=False, apply_mult=True): structure = self.data_dimensions.structure package_dim = self.data_dimensions.package_dim for cellid in model_grid.get_all_model_cells(): + first_item = self.layer_storage.first_item() data_line = (cellid,) + ( - self.layer_storage.first_item().data_const_value, + first_item.data_const_value, ) if len(structure.data_item_structures) > 2: # append None any expected optional data @@ -824,9 +826,8 @@ def append_data(self, data): np.rec.array(internal_data_list, self._recarray_type_list) ) else: - if len(self.layer_storage.first_item().internal_data[0]) > len( - data[0] - ): + first_item = self.layer_storage.first_item() + if len(first_item.internal_data[0]) > len(data[0]): # Add placeholders to data self._add_placeholders(data) self.set_data( @@ -1249,7 +1250,7 @@ def make_tuple_cellids(self, data): new_data = [] current_cellid = () - for index, line in enumerate(data): + for line in data: new_line = [] for item, is_cellid in zip(line, self.recarray_cellid_list_ex): if is_cellid: @@ -1434,7 +1435,9 @@ def store_external( self._data_path, self._stress_period, ) - file_access.write_text_file(data, fp, data_type, data_size) + file_access.write_text_file( + data, fp, data_type, data_size, + ) self.layer_storage[layer_new].factor = multiplier self.layer_storage[layer_new].internal_data = None self.layer_storage[layer_new].data_const_value = None @@ -1663,7 +1666,12 @@ def internal_to_external( ) def resolve_shape_list( - self, data_item, repeat_count, current_key, data_line, cellid_size=None + self, + data_item, + repeat_count, + current_key, + data_line, + cellid_size=None, ): struct = self.data_dimensions.structure try: @@ -1742,18 +1750,15 @@ def process_internal_line(self, arr_line): index = 1 while index < len(arr_line): if isinstance(arr_line[index], str): - if arr_line[index].lower() == "factor" and index + 1 < len( - arr_line - ): + word = arr_line[index].lower() + if word == "factor" and index + 1 < len(arr_line): multiplier = convert_data( arr_line[index + 1], self.data_dimensions, self._data_type, ) index += 2 - elif arr_line[index].lower() == "iprn" and index + 1 < len( - arr_line - ): + elif word == "iprn" and index + 1 < len(arr_line): print_format = arr_line[index + 1] index += 2 else: @@ -1808,9 +1813,8 @@ def process_open_close_line(self, arr_line, layer, store=True): ) while index < len(arr_line): if isinstance(arr_line[index], str): - if arr_line[index].lower() == "factor" and index + 1 < len( - arr_line - ): + word = arr_line[index].lower() + if word == "factor" and index + 1 < len(arr_line): try: multiplier = convert_data( arr_line[index + 1], @@ -1840,20 +1844,13 @@ def process_open_close_line(self, arr_line, layer, store=True): ex, ) index += 2 - elif arr_line[index].lower() == "iprn" and index + 1 < len( - arr_line - ): + elif word == "iprn" and index + 1 < len(arr_line): print_format = arr_line[index + 1] index += 2 - elif arr_line[index].lower() == "data" and index + 1 < len( - arr_line - ): + elif word == "data" and index + 1 < len(arr_line): data = arr_line[index + 1] index += 2 - elif ( - arr_line[index].lower() == "binary" - or arr_line[index].lower() == "(binary)" - ): + elif word == "binary" or word == "(binary)": binary = True index += 1 else: @@ -1966,16 +1963,16 @@ def _verify_list(self, data): for index in range( 0, min(data_line_len, len(self._recarray_type_list)) ): + datadim = self.data_dimensions if ( self._recarray_type_list[index][0] == "cellid" - and self.data_dimensions.get_model_dim(None).model_name - is not None + and datadim.get_model_dim(None).model_name is not None and data_line[index] is not None ): # this is a cell id. verify that it contains the # correct number of integers if cellid_size is None: - model_grid = self.data_dimensions.get_model_grid() + model_grid = datadim.get_model_grid() cellid_size = ( model_grid.get_num_spatial_coordinates() ) @@ -2423,12 +2420,12 @@ def build_type_list( and data is not None ): idx = 1 - data_line_max_size = self._get_max_data_line_size(data) + line_max_size = self._get_max_data_line_size(data) type_list = self.resolve_typelist(data) - while len(type_list) < data_line_max_size: - # keystrings at the end of a line can contain items - # of variable length. assume everything at the - # end of the data line is related to the last + while len(type_list) < line_max_size: + # keystrings at the end of a line can contain + # items of variable length. assume everything at + # the end of the data line is related to the last # keystring name = "{}_{}".format(ks_data_item.name, idx) self._append_type_lists( @@ -2468,14 +2465,14 @@ def build_type_list( model_dim = self.data_dimensions.get_model_dim( None ) - expression_array = model_dim.build_shape_expression( + exp_array = model_dim.build_shape_expression( data_item.shape ) if ( - isinstance(expression_array, list) - and len(expression_array) == 1 + isinstance(exp_array, list) + and len(exp_array) == 1 ): - exp = expression_array[0] + exp = exp_array[0] resolved_shape = [ model_dim.resolve_exp(exp, nseg) ] @@ -2529,8 +2526,8 @@ def build_type_list( # A cellid is a single entry (tuple) in the # recarray. Adjust dimensions accordingly. data_dim = self.data_dimensions - model_grid = data_dim.get_model_grid() - size = model_grid.get_num_spatial_coordinates() + grid = data_dim.get_model_grid() + size = grid.get_num_spatial_coordinates() data_item.remove_cellid(resolved_shape, size) for index in range(0, resolved_shape[0]): if resolved_shape[0] > 1: @@ -2562,22 +2559,22 @@ def _append_type_lists(self, name, data_type, iscellid): cellid_size = model_grid.get_num_spatial_coordinates() # determine header for different grid types if cellid_size == 1: - self._recarray_type_list_ex.append((name, int)) + self._do_ex_list_append(name, int, iscellid) elif cellid_size == 2: - self._recarray_type_list_ex.append(("layer", int)) - self._recarray_type_list_ex.append(("cell2d_num", int)) + self._do_ex_list_append("layer", int, iscellid) + self._do_ex_list_append("cell2d_num", int, iscellid) else: - self._recarray_type_list_ex.append(("layer", int)) - self._recarray_type_list_ex.append(("row", int)) - self._recarray_type_list_ex.append(("column", int)) - - for index in range(0, cellid_size): - self.recarray_cellid_list_ex.append(iscellid) + self._do_ex_list_append("layer", int, iscellid) + self._do_ex_list_append("row", int, iscellid) + self._do_ex_list_append("column", int, iscellid) else: - self._recarray_type_list_ex.append((name, data_type)) - self.recarray_cellid_list_ex.append(iscellid) + self._do_ex_list_append(name, data_type, iscellid) return iscellid + def _do_ex_list_append(self, name, data_type, iscellid): + self._recarray_type_list_ex.append((name, data_type)) + self.recarray_cellid_list_ex.append(iscellid) + @staticmethod def _calc_data_size(data, count_to=None, current_length=None): if current_length is None: diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 0077996871..01d7f97890 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -733,11 +733,11 @@ def dtype( model.simulation_data.debug, ) if aux_vars is not None: - if len(aux_vars) > 0 and ( - isinstance(aux_vars[0], list) - or isinstance(aux_vars[0], tuple) - ): - aux_vars = aux_vars[0] + if len(aux_vars) > 0: + if isinstance(aux_vars[0], list) or isinstance( + aux_vars[0], tuple + ): + aux_vars = aux_vars[0] for aux_var in aux_vars: type_list.append((aux_var, object)) if boundnames: @@ -768,7 +768,7 @@ def empty( ) # get data storage - data_struct, data_dimensions = self._get_data_dimensions(model) + data_struct = self._get_data_dimensions(model)[0] data_type = data_struct.get_datatype() # build recarray From db0937006303832ad6ec304a32a22acc6097687e Mon Sep 17 00:00:00 2001 From: spaulins-usgs Date: Thu, 27 Aug 2020 12:04:00 -0700 Subject: [PATCH 6/6] fix(formatting): using the latest version of black --- flopy/mf6/data/mfdatalist.py | 5 ++++- flopy/mf6/data/mfdatastorage.py | 12 ++++++------ flopy/mf6/data/mfdatautil.py | 7 ++++++- 3 files changed, 16 insertions(+), 8 deletions(-) diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index e4ce725c6d..88027883ca 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -206,7 +206,10 @@ def to_array(self, kper=0, mask=False): model_grid.num_columns(), ) elif model_grid._grid_type.value == 2: - shape = (model_grid.num_layers(), model_grid.num_cells_per_layer()) + shape = ( + model_grid.num_layers(), + model_grid.num_cells_per_layer(), + ) else: shape = (model_grid.num_cells_per_layer(),) diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py index 0219a96637..88fcd2a752 100644 --- a/flopy/mf6/data/mfdatastorage.py +++ b/flopy/mf6/data/mfdatastorage.py @@ -1436,7 +1436,10 @@ def store_external( self._stress_period, ) file_access.write_text_file( - data, fp, data_type, data_size, + data, + fp, + data_type, + data_size, ) self.layer_storage[layer_new].factor = multiplier self.layer_storage[layer_new].internal_data = None @@ -2465,11 +2468,8 @@ def build_type_list( model_dim = self.data_dimensions.get_model_dim( None ) - - expression_array = ( - model_dim.build_shape_expression( - data_item.shape - ) + exp_array = model_dim.build_shape_expression( + data_item.shape ) if ( isinstance(exp_array, list) diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 01d7f97890..dcc30cc01c 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -764,7 +764,12 @@ def empty( # get type list type_list = self.dtype( - model, aux_vars, boundnames, nseg, timeseries, cellid_expanded, + model, + aux_vars, + boundnames, + nseg, + timeseries, + cellid_expanded, ) # get data storage