diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 142fed140d..f5eba2f233 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -1450,3 +1450,64 @@ def test_geo_dataframe(structured_grid, vertex_grid, unstructured_grid): raise AssertionError( f"Cell vertices incorrect for node={node}" ) + + +def test_unstructured_iverts_cleanup(): + grid = GridCases.structured_small() + + # begin building unstructured grid information + top = grid.top.ravel() + botm = grid.botm[0].ravel() + idomain = np.ones(botm.shape, dtype=int) + + # build iac and ja + neighbors = grid.neighbors(method="rook", reset=True) + iac, ja = [], [] + for cell, neigh in neighbors.items(): + iac.append(len(neigh) + 1) + ja.extend( + [ + cell, + ] + + neigh + ) + + # build iverts and verts without using shared vertices + verts, iverts = [], [] + xverts, yverts = grid.cross_section_vertices + ivt = 0 + for cid, xvs in enumerate(xverts): + yvs = yverts[cid] + ivts = [] + for ix, vert in enumerate(xvs[:-1]): + ivts.append(ivt) + verts.append([ivt, vert, yvs[ix]]) + ivt += 1 + + ivts.append(ivts[0]) + iverts.append(ivts) + + ugrid = UnstructuredGrid( + vertices=verts, + iverts=iverts, + xcenters=grid.xcellcenters.ravel(), + ycenters=grid.ycellcenters.ravel(), + iac=iac, + ja=ja, + top=top, + botm=botm, + idomain=idomain, + ) + + if ugrid.nvert != (grid.ncpl * 4): + raise AssertionError( + "UnstructuredGrid is being built incorrectly for test case" + ) + + cleaned_vert_num = (grid.nrow + 1) * (grid.ncol + 1) + clean_ugrid = ugrid.clean_iverts() + + if clean_ugrid.nvert != cleaned_vert_num: + raise AssertionError( + "Improper number of vertices for cleaned 'shared' iverts" + ) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index c7756589d1..c7b3c4b53e 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -26,7 +26,7 @@ def test_structured_model_splitter(function_tmpdir): for row in range(modelgrid.nrow): if row != 0 and row % 2 == 0: ncol += 1 - array[row, ncol:] = 2 + array[row, ncol:] = 100 mfsplit = Mf6Splitter(sim) new_sim = mfsplit.split_model(array) @@ -37,13 +37,13 @@ def test_structured_model_splitter(function_tmpdir): original_heads = gwf.output.head().get_alldata()[-1] - ml0 = new_sim.get_model("freyberg_1") - ml1 = new_sim.get_model("freyberg_2") + ml0 = new_sim.get_model("freyberg_001") + ml1 = new_sim.get_model("freyberg_100") heads0 = ml0.output.head().get_alldata()[-1] heads1 = ml1.output.head().get_alldata()[-1] - new_heads = mfsplit.reconstruct_array({1: heads0, 2: heads1}) + new_heads = mfsplit.reconstruct_array({1: heads0, 100: heads1}) err_msg = "Heads from original and split models do not match" np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg) @@ -691,3 +691,171 @@ def test_idomain_none(function_tmpdir): err_msg = "Heads from original and split models do not match" np.testing.assert_allclose(new_head, head, atol=1e-07, err_msg=err_msg) + + +@requires_exe("mf6") +def test_unstructured_complex_disu(function_tmpdir): + sim_path = function_tmpdir + split_sim_path = sim_path / "model_split" + + # build the simulation structure + sim = flopy.mf6.MFSimulation(sim_ws=sim_path) + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + mname = "disu_model" + gwf = flopy.mf6.ModflowGwf(sim, modelname=mname) + + # start structured and then create a USG from it + nlay = 1 + nrow = 10 + ncol = 10 + delc = np.ones((nrow,)) + delr = np.ones((ncol,)) + top = np.ones((nrow, ncol)) + botm = np.zeros((nlay, nrow, ncol)) + idomain = np.ones(botm.shape, dtype=int) + idomain[0, 1, 4] = 0 + idomain[0, 8, 5] = 0 + + grid = flopy.discretization.StructuredGrid( + delc=delc, delr=delr, top=top, botm=botm, idomain=idomain + ) + + # build the USG connection information + neighbors = grid.neighbors(method="rook", reset=True) + iac, ja, ihc, cl12, hwva, angldegx = [], [], [], [], [], [] + for cell, neigh in neighbors.items(): + iac.append(len(neigh) + 1) + ihc.extend( + [ + 1, + ] + * (len(neigh) + 1) + ) + ja.extend( + [ + cell, + ] + + neigh + ) + cl12.extend( + [ + 0, + ] + + [ + 1, + ] + * len(neigh) + ) + hwva.extend( + [ + 0, + ] + + [ + 1, + ] + * len(neigh) + ) + adx = [ + 0, + ] + for n in neigh: + ev = cell - n + if ev == -1 * ncol: + adx.append(270) + elif ev == ncol: + adx.append(90) + elif ev == -1: + adx.append(0) + else: + adx.append(180) + angldegx.extend(adx) + + # build iverts and verts. Do not use shared iverts and mess with verts a + # tiny bit + verts, cell2d = [], [] + xverts, yverts = grid.cross_section_vertices + xcenters = grid.xcellcenters.ravel() + ycenters = grid.ycellcenters.ravel() + ivert = 0 + for cell_num, xvs in enumerate(xverts): + if (cell_num - 3) % 10 == 0: + xvs[2] -= 0.001 + xvs[3] -= 0.001 + yvs = yverts[cell_num] + + c2drec = [cell_num, xcenters[cell_num], ycenters[cell_num], len(xvs)] + for ix, vert in enumerate(xvs[:-1]): + c2drec.append(ivert) + verts.append([ivert, vert, yvs[ix]]) + ivert += 1 + + c2drec.append(c2drec[4]) + cell2d.append(c2drec) + + nodes = len(cell2d) + nja = len(ja) + nvert = len(verts) + + dis = flopy.mf6.ModflowGwfdisu( + gwf, + nodes=nodes, + nja=nja, + nvert=nvert, + top=np.ravel(grid.top), + bot=np.ravel(grid.botm), + area=np.ones((nodes,)), + idomain=grid.idomain.ravel(), + iac=iac, + ja=ja, + ihc=ihc, + cl12=cl12, + hwva=hwva, + angldegx=angldegx, + vertices=verts, + cell2d=cell2d, + ) + + # build npf, ic, CHD, OC package + npf = flopy.mf6.ModflowGwfnpf(gwf) + ic = flopy.mf6.ModflowGwfic(gwf) + + spd = [] + for i in range(nrow): + spd.append((0 + (i * 10), 0.9)) + spd.append((9 + (i * 10), 0.5)) + + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=spd) + + spd = {0: [("HEAD", "LAST")]} + oc = flopy.mf6.ModflowGwfoc( + gwf, head_filerecord=f"{mname}.hds", saverecord=spd + ) + + sim.write_simulation() + sim.run_simulation() + + heads = gwf.output.head().get_alldata()[-1] + + array = np.zeros((nrow, ncol), dtype=int) + array[:, 5:] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array) + + new_sim.set_sim_path(split_sim_path) + new_sim.write_simulation() + new_sim.run_simulation() + + gwf0 = new_sim.get_model(f"{mname}_0") + gwf1 = new_sim.get_model(f"{mname}_1") + + heads0 = gwf0.output.head().get_alldata()[-1] + heads1 = gwf1.output.head().get_alldata()[-1] + + new_heads = mfsplit.reconstruct_array({0: heads0, 1: heads1}) + + diff = np.abs(heads - new_heads) + if np.max(diff) > 1e-07: + raise AssertionError("Reconstructed head results outside of tolerance") diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index 24d6b6cf73..ed0a201930 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -131,6 +131,7 @@ def __init__( angrot=0.0, iac=None, ja=None, + cell2d=None, **kwargs, ): super().__init__( @@ -146,6 +147,11 @@ def __init__( angrot=angrot, **kwargs, ) + if cell2d is not None: + # modflow 6 DISU + xcenters = np.array([i[1] for i in cell2d]) + ycenters = np.array([i[2] for i in cell2d]) + iverts = [list(t)[4:] for t in cell2d] # if any of these are None, then the grid is not valid self._vertices = vertices @@ -239,6 +245,31 @@ def nnodes(self): def nvert(self): return len(self._vertices) + @property + def cell2d(self): + if self.is_valid: + ncenters = len(self._xc) + is_layered = False + if ncenters != self.nnodes and ncenters / self.nnodes % 0: + is_layered = True + + ix_adj = 0 + cell2d = [] + for ix in range(self.nnodes): + if is_layered: + if ix % self.ncpl[0] == 0 and ix != 0: + ix_adj += self.ncpl[0] + iverts = self._iverts[ix - ix_adj] + c2drec = [ + ix, + self._xc[ix - ix_adj], + self._yc[ix - ix_adj], + len(iverts), + ] + c2drec.extend(iverts) + cell2d.append(c2drec) + return cell2d + @property def iverts(self): if self._iverts is not None: @@ -566,6 +597,45 @@ def geo_dataframe(self): gdf = super().geo_dataframe(polys) return gdf + def neighbors(self, node=None, **kwargs): + """ + Method to get nearest neighbors of a cell + + Parameters + ---------- + node : int + model grid node number + + ** kwargs: + method : str + "iac" for specified connections from the DISU package + "rook" for shared edge neighbors + "queen" for shared vertex neighbors + reset : bool + flag to reset the neighbor calculation + + Returns + ------- + list or dict : list of cell node numbers or dict of all cells and + neighbors + """ + method = kwargs.pop("method", None) + reset = kwargs.pop("reset", False) + if method == "iac": + if self._neighbors is None or reset: + neighors = {} + idx0 = 0 + for node, ia in enumerate(self._iac): + idx1 = idx0 + ia + neighors[node] = list(self._ja[idx0 + 1 : idx1]) + self._neighbors = neighors + if node is not None: + return self._neighbors[node] + else: + return self._neighbors + else: + return super().neighbors(node=node, method=method, reset=reset) + def convert_grid(self, factor): """ Method to scale the model grid based on user supplied scale factors @@ -599,6 +669,66 @@ def convert_grid(self, factor): "Grid is not complete and cannot be converted" ) + def clean_iverts(self, inplace=False): + """ + Method to clean up duplicated iverts/verts when vertex information + is supplied in the unstructured grid. + + Parameters: + ---------- + inplace : bool + flag to clean and reset iverts in the current modelgrid object. + Default is False and returns a new modelgrid object + + Returns + ------- + UnstructuredGrid or None + """ + if self.is_valid: + vset = {} + for rec in self._vertices: + vert = (rec[1], rec[2]) + if vert in vset: + vset[vert].add(rec[0]) + else: + vset[vert] = { + rec[0], + } + + cnt = 0 + ivert_remap = {} + vertices = [] + for (xv, yv), iverts in vset.items(): + for iv in iverts: + ivert_remap[iv] = cnt + vertices.append((cnt, xv, yv)) + cnt += 1 + + iverts = [[ivert_remap[v] for v in ivs] for ivs in self.iverts] + if inplace: + self._vertices = vertices + self._iverts = iverts + self._require_cache_updates() + else: + return UnstructuredGrid( + vertices, + iverts=iverts, + xcenters=self._xc, + ycenters=self._yc, + top=self._top, + botm=self._botm, + idomain=self._idomain, + lenuni=self.lenuni, + ncpl=self._ncpl, + crs=self._crs, + prjfile=self._prjfile, + xoff=self.xoffset, + yoff=self.yoffset, + angrot=self.angrot, + iac=self._iac, + ja=self._ja, + ) + def intersect(self, x, y, z=None, local=False, forgive=False): """ Get the CELL2D number of a point with coordinates x and y diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index 70c62e28e4..cd12f6ca0a 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -445,7 +445,9 @@ def modelgrid(self): if ncpl is None: ncpl = np.array([dis.nodes.get_data()], dtype=int) cell2d = dis.cell2d.array - idomain = np.ones(dis.nodes.array, np.int32) + idomain = dis.idomain.array + if idomain is None: + idomain = np.ones(dis.nodes.array, dtype=int) if cell2d is None: if ( self.simulation.simulation_data.verbosity_level.value @@ -455,13 +457,7 @@ def modelgrid(self): "WARNING: cell2d information missing. Functionality of " "the UnstructuredGrid will be limited." ) - iverts = None - xcenters = None - ycenters = None - else: - iverts = [list(i)[4:] for i in cell2d] - xcenters = dis.cell2d.array["xc"] - ycenters = dis.cell2d.array["yc"] + vertices = dis.vertices.array if vertices is None: if ( @@ -478,9 +474,7 @@ def modelgrid(self): self._modelgrid = UnstructuredGrid( vertices=vertices, - iverts=iverts, - xcenters=xcenters, - ycenters=ycenters, + cell2d=cell2d, top=dis.top.array, botm=dis.bot.array, idomain=idomain, diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index b16ba93b6b..a73daf27af 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -152,6 +152,7 @@ def __init__(self, sim, modelname=None): self._connection = None self._uconnection = None self._usg_metadata = None + self._has_angldegx = False self._connection_ivert = None self._model_dict = None self._ivert_vert_remap = None @@ -170,6 +171,8 @@ def __init__(self, sim, modelname=None): self._maw_remaps = {} self._allow_splitting = True + self._fdigits = 1 + @property def new_simulation(self): """ @@ -216,6 +219,7 @@ def switch_models(self, modelname, remap_nodes=False): self._connection = None self._uconnection = None self._usg_metadata = None + self._has_angldegx = False self._connection_ivert = None self._ivert_vert_remap = None self._sfr_mover_connections = [] @@ -657,8 +661,7 @@ def _remap_nodes(self, array): f"{bad_keys} are not in the active model extent; " f"please adjust the model splitting array" ) - - if self._modelgrid.iverts is None: + if self._modelgrid.grid_type == "unstructured": self._map_iac_ja_connections() else: self._connection = self._modelgrid.neighbors( @@ -703,14 +706,7 @@ def _remap_nodes(self, array): for m in np.unique(array): cells = np.asarray(array == m).nonzero()[0] - mapping = np.zeros( - ( - len( - cells, - ) - ), - dtype=int, - ) + mapping = np.zeros((len(cells),), dtype=int) mapping[:] = cells grid_info[m] = [(len(cells),), None, None, mapping] @@ -760,7 +756,7 @@ def _remap_nodes(self, array): if cmdl == mdl: if nnode in new_connections[mdl]["internal"]: new_connections[mdl]["internal"][nnode].append(cnnode) - if self._connection_ivert is None: + if self._uconnection is not None: usg_meta[mdl][nnode]["ihc"].append( int(self._uconnection[node]["ihc"][ix + 1]) ) @@ -770,10 +766,14 @@ def _remap_nodes(self, array): usg_meta[mdl][nnode]["hwva"].append( self._uconnection[node]["hwva"][ix + 1] ) + if self._has_angldegx: + usg_meta[mdl][nnode]["angldegx"].append( + self._uconnection[node]["angldegx"][ix + 1] + ) else: new_connections[mdl]["internal"][nnode] = [cnnode] - if self._connection_ivert is None: + if self._uconnection is not None: usg_meta[mdl][nnode] = { "ihc": [ self._uconnection[node]["ihc"][0], @@ -788,14 +788,20 @@ def _remap_nodes(self, array): self._uconnection[node]["hwva"][ix + 1], ], } + if self._has_angldegx: + usg_meta[mdl][nnode]["angldegx"] = [ + self._uconnection[node]["angldegx"][0], + self._uconnection[node]["angldegx"][ + ix + 1 + ], + ] else: if nnode in new_connections[mdl]["external"]: new_connections[mdl]["external"][nnode].append( (cmdl, cnnode) ) - if self._connection_ivert is not None: - tmp = self._connection_ivert[node] + if self._uconnection is None: exchange_meta[mdl][nnode][cnnode] = [ node, cnode, @@ -813,7 +819,7 @@ def _remap_nodes(self, array): new_connections[mdl]["external"][nnode] = [ (cmdl, cnnode) ] - if self._connection_ivert is not None: + if self._uconnection is None: exchange_meta[mdl][nnode] = { cnnode: [ node, @@ -832,7 +838,7 @@ def _remap_nodes(self, array): ] } - if self._modelgrid.grid_type == "vertex": + if self._modelgrid.grid_type in ("vertex", "unstructured"): self._map_verts_iverts(array) self._new_connections = new_connections @@ -854,16 +860,23 @@ def _map_iac_ja_connections(self): cl12 = self._model.disu.cl12.array ihc = self._model.disu.ihc.array hwva = self._model.disu.hwva.array + angldegx = self._model.disu.angldegx.array + if angldegx is not None: + self._has_angldegx = True idx0 = 0 for ia in iac: idx1 = idx0 + ia cn = ja[idx0 + 1 : idx1] - conn[ja[idx0]] = cn + conn[ja[idx0]] = list(cn) uconn[ja[idx0]] = { - "cl12": cl12[idx0:idx1], - "ihc": ihc[idx0:idx1], - "hwva": hwva[idx0:idx1], + "cl12": list(cl12[idx0:idx1]), + "ihc": list(ihc[idx0:idx1]), + "hwva": list(hwva[idx0:idx1]), } + + if self._has_angldegx: + uconn[ja[idx0]]["angldegx"] = list(angldegx[idx0:idx1]) + idx0 = idx1 self._connection = conn @@ -881,6 +894,8 @@ def _map_verts_iverts(self, array): """ iverts = self._modelgrid.iverts verts = self._modelgrid.verts + if iverts is None: + return ivlut = {mkey: {} for mkey in np.unique(array)} for mkey in np.unique(array): @@ -1013,9 +1028,7 @@ def _remap_filerecords(self, item, value, mapped_data): for mdl in mapped_data.keys(): if mapped_data[mdl]: new_val = value.split(".") - new_val = ( - f"{'.'.join(new_val[0:-1])}_{mdl}.{new_val[-1]}" - ) + new_val = f"{'.'.join(new_val[0:-1])}_{mdl :0{self._fdigits}d}.{new_val[-1]}" mapped_data[mdl][item] = new_val return mapped_data @@ -1034,7 +1047,7 @@ def _remap_disu(self, mapped_data): """ for mkey, metadata in self._usg_metadata.items(): - iac, ja, ihc, cl12, hwva = [], [], [], [], [] + iac, ja, ihc, cl12, hwva, angldegx = [], [], [], [], [], [] for node, params in metadata.items(): conns = [node] + self._new_connections[mkey]["internal"][node] iac.append(len(conns)) @@ -1042,6 +1055,8 @@ def _remap_disu(self, mapped_data): ihc.extend(params["ihc"]) cl12.extend(params["cl12"]) hwva.extend(params["hwva"]) + if self._has_angldegx: + angldegx.extend(params["angldegx"]) assert np.sum(iac) == len(ja) @@ -1051,6 +1066,8 @@ def _remap_disu(self, mapped_data): mapped_data[mkey]["ihc"] = ihc mapped_data[mkey]["cl12"] = cl12 mapped_data[mkey]["hwva"] = hwva + if self._has_angldegx: + mapped_data[mkey]["angldegx"] = angldegx return mapped_data @@ -1195,7 +1212,7 @@ def _remap_array(self, item, mfarray, mapped_data, **kwargs): else: # external array tmp = fnames[lay].split(".") - filename = f"{'.'.join(tmp[:-1])}.{mkey}.{tmp[-1]}" + filename = f"{'.'.join(tmp[:-1])}.{mkey :0{self._fdigits}d}.{tmp[-1]}" cr = { "filename": filename, @@ -1280,7 +1297,7 @@ def _remap_mflist( if how == 3 and new_recarray is not None: tmp = fname.split(".") - filename = f"{'.'.join(tmp[:-1])}.{mkey}.{tmp[-1]}" + filename = f"{'.'.join(tmp[:-1])}.{mkey :0{self._fdigits}d}.{tmp[-1]}" new_recarray = { "data": new_recarray, @@ -1329,7 +1346,7 @@ def _remap_adv_transport(self, package, item, pkg_remap, mapped_data): ) spd[per] = new_recarray - flow_package_const[-2] += f"_{mkey}" + flow_package_const[-2] += f"_{mkey :0{self._fdigits}d}" new_flow_package_name = ".".join(flow_package_const) mapped_data[mkey]["packagedata"] = new_packagedata mapped_data[mkey]["perioddata"] = spd @@ -2218,7 +2235,7 @@ def _remap_fmi(self, package, mapped_data): new_fnames = [] for fname in fnames: new_val = fname.split(".") - new_val = f"{'.'.join(new_val[0:-1])}_{mkey}.{new_val[-1]}" + new_val = f"{'.'.join(new_val[0:-1])}_{mkey :0{self._fdigits}d}.{new_val[-1]}" new_fnames.append(new_val) new_packagedata = packagedata.copy() @@ -2562,11 +2579,11 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): if isinstance(ofile, (list, tuple)): fname = ofile[0] tmp = fname.split(".") - tmp[-2] += f"_{mkey}" + tmp[-2] += f"_{mkey :0{self._fdigits}d}" ofile[0] = ".".join(tmp) else: tmp = ofile.split(".") - tmp[-2] += f"_{mkey}" + tmp[-2] += f"_{mkey :0{self._fdigits}d}" ofile = ".".join(tmp) if pkg_type is None: @@ -2804,6 +2821,9 @@ def _remap_package(self, package, ismvr=False): ), ): for item, value in package.__dict__.items(): + if item in ("nja", "ja", "cl12", "ihc", "hwva", "angldegx"): + continue + if item in ("delr", "delc"): for mkey, d in self._grid_info.items(): if item == "delr": @@ -2832,13 +2852,27 @@ def _remap_package(self, package, ismvr=False): elif item == "iac": mapped_data = self._remap_disu(mapped_data) - break elif item == "xorigin": for mkey in self._model_dict.keys(): for k, v in self._offsets[mkey].items(): mapped_data[mkey][k] = v + elif item in ("vertices", "cell2d"): + if value.array is not None: + if item == "cell2d": + mapped_data = self._remap_cell2d( + item, value, mapped_data + ) + else: + for mkey in self._model_dict.keys(): + mapped_data[mkey][item] = ( + self._ivert_vert_remap[mkey][item] + ) + mapped_data[mkey]["nvert"] = len( + self._ivert_vert_remap[mkey][item] + ) + elif isinstance(value, mfdataarray.MFArray): mapped_data = self._remap_array(item, value, mapped_data) @@ -3066,19 +3100,26 @@ def _create_exchanges(self): if self._modelgrid.grid_type == "unstructured": # use existing connection information aux = False + grid_dict = {i: m.modelgrid for i, m in self._model_dict.items()} for m0, model in self._model_dict.items(): + grid0 = grid_dict[m0] exg_nodes = self._new_connections[m0]["external"] for m1 in nmodels: + grid1 = grid_dict[m1] if m1 in built: continue if m1 == m0: continue exchange_data = [] for node0, exg_list in exg_nodes.items(): + if grid0.idomain[node0] < 1: + continue for exg in exg_list: if exg[0] != m1: continue node1 = exg[-1] + if grid1.idomain[node1] < 1: + continue exg_meta0 = self._exchange_metadata[m0][node0][ node1 ] @@ -3104,7 +3145,7 @@ def _create_exchanges(self): exgmnameb=mname1, nexg=len(exchange_data), exchangedata=exchange_data, - filename=f"sim_{m0}_{m1}.{extension}", + filename=f"sim_{m0 :0{self._fdigits}d}_{m1 :0{self._fdigits}d}.{extension}", newton=newton, xt3d=xt3d, ) @@ -3242,7 +3283,7 @@ def _create_exchanges(self): auxiliary=["ANGLDEGX", "CDIST"], nexg=len(exchange_data), exchangedata=exchange_data, - filename=f"sim_{m0}_{m1}.{extension}", + filename=f"sim_{m0 :0{self._fdigits}d}_{m1 :0{self._fdigits}d}.{extension}", newton=newton, xt3d=xt3d, ) @@ -3286,6 +3327,11 @@ def split_model(self, array): "is part of a split simulation" ) + # set number formatting string for file paths + array = np.array(array).astype(int) + s = str(np.max(array)) + self._fdigits = len(s) + self._remap_nodes(array) if self._new_sim is None: @@ -3306,7 +3352,7 @@ def split_model(self, array): mdl_cls = PackageContainer.model_factory(self._model_type) self._model_dict[mkey] = mdl_cls( self._new_sim, - modelname=f"{self._modelname}_{mkey}", + modelname=f"{self._modelname}_{mkey :0{self._fdigits}d}", **nam_options, ) @@ -3320,9 +3366,3 @@ def split_model(self, array): epaks = self._create_exchanges() return self._new_sim - - -# todo: development notes: -# Then set up checks for model splitting -# (ex. doesn't parallel a fault, doesn't cut through a lake, -# active cells in modelgrid...) diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py index fe91551213..4fabf0482b 100644 --- a/flopy/utils/rasters.py +++ b/flopy/utils/rasters.py @@ -537,6 +537,13 @@ def resample_to_grid( arr = self.get_array(band, masked=True) arr = arr.flatten() + # filter out nan values from the original dataset + if np.isnan(np.sum(arr)): + idx = np.isfinite(arr) + rxc = rxc[idx] + ryc = ryc[idx] + arr = arr[idx] + # step 3: use griddata interpolation to snap to grid data = griddata( (rxc, ryc),