diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 2cd4e9096..8861789fa 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -285,6 +285,28 @@ def test_unstructured_grid_get_cell_vertices(): assert v1 == v4, "Positional and kwarg should match" +def test_get_cell_iverts(): + nrow, ncol = 4, 5 + grid = StructuredGrid(delr=np.ones(ncol), delc=np.ones(nrow)) + + iverts = grid.get_cell_iverts() + assert iverts.shape == (grid.ncpl, 4) + assert np.array_equal(iverts.tolist(), grid.iverts) + + nodes = [0, 7, grid.ncpl - 1] + assert np.array_equal(grid.get_cell_iverts(nodes), iverts[nodes]) + assert np.array_equal(grid.get_cell_iverts(np.array(nodes)), iverts[nodes]) + assert np.array_equal( + grid.get_cell_iverts(np.array(nodes).reshape((-1, 1))), iverts[nodes] + ) + assert np.array_equal(grid.get_cell_iverts(3), iverts[[3]]) + + # vertices are ordered clockwise from the upper left corner of a cell + for node, iv in enumerate(iverts): + i, j = divmod(node, ncol) + assert list(iv) == grid._build_structured_iverts(i, j) + + def test_get_lrc_get_node(): nlay, nrow, ncol = 3, 4, 5 nnodes = nlay * nrow * ncol diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 60bbb8330..14556d44a 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -138,6 +138,321 @@ def test_unstructured_model_splitter(function_tmpdir): np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg) +@requires_exe("mf6") +def test_structured_to_disv_model_splitter(function_tmpdir): + sim_path = get_example_data_path() / "mf6-freyberg" + split_path = function_tmpdir / "split_model" + + sim = MFSimulation.load(sim_ws=sim_path) + sim.set_sim_path(function_tmpdir) + sim.write_simulation() + sim.run_simulation() + + gwf = sim.get_model() + modelgrid = gwf.modelgrid + idomain = modelgrid.idomain.reshape((modelgrid.nlay, modelgrid.ncpl)) + + array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int) + array[modelgrid.nrow // 2 :, :] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array, to_disv=True) + + new_sim.set_sim_path(split_path) + new_sim.write_simulation() + new_sim.run_simulation() + + ml0 = new_sim.get_model("freyberg_0") + ml1 = new_sim.get_model("freyberg_1") + + for ml in (ml0, ml1): + if ml.modelgrid.grid_type != "vertex": + raise AssertionError("Split model is not a DISV model") + + # every stack of cells that is active in a layer, and no other, is split + nactive = np.count_nonzero(np.any(idomain != 0, axis=0)) + if ml0.modelgrid.ncpl + ml1.modelgrid.ncpl != nactive: + raise AssertionError("Split models do not include all active cells") + + original_heads = gwf.output.head().get_alldata()[-1] + heads0 = ml0.output.head().get_alldata()[-1] + heads1 = ml1.output.head().get_alldata()[-1] + + new_heads = mfsplit.reconstruct_array({0: heads0, 1: heads1}) + + idx = modelgrid.idomain != 0 + err_msg = "Heads from original and split models do not match" + np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg) + + +@requires_exe("mf6") +def test_structured_to_disv_idomain(function_tmpdir): + nlay, nrow, ncol = 3, 10, 10 + idomain = np.ones((nlay, nrow, ncol), dtype=int) + idomain[:, 0:2, 0:2] = 0 # inactive in every layer, excluded from the split + idomain[0, 5:, 5:] = 0 # inactive in one layer, carried into the split + idomain[2, 0:3, 7:] = 0 + idomain[1, 3:7, 3:7] = -1 # vertical passthrough, carried into the split + + sim = flopy.mf6.MFSimulation(sim_name="ml", sim_ws=function_tmpdir, exe_name="mf6") + flopy.mf6.ModflowTdis(sim) + flopy.mf6.ModflowIms(sim, complexity="simple") + gwf = flopy.mf6.ModflowGwf(sim, modelname="ml", save_flows=True) + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=100.0, + delc=100.0, + top=30.0, + botm=[20.0, 10.0, 0.0], + idomain=idomain, + ) + flopy.mf6.ModflowGwfnpf(gwf, k=1.0) + flopy.mf6.ModflowGwfic(gwf, strt=25.0) + chd = [[(0, i, 2), 25.0] for i in range(2, nrow)] + chd += [[(0, i, ncol - 1), 20.0] for i in range(0, 5)] + flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd) + flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord="ml.hds", + saverecord=[("HEAD", "ALL")], + ) + sim.write_simulation() + sim.run_simulation() + + original_heads = gwf.output.head().get_alldata()[-1] + + # split on a raveled array that covers every cell of the DIS grid + array = np.zeros((nrow, ncol), dtype=int) + array[:, ncol // 2 :] = 1 + array = array.ravel() + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array, to_disv=True) + new_sim.set_sim_path(function_tmpdir / "split_model") + new_sim.write_simulation() + new_sim.run_simulation() + + heads = {} + ncpl = 0 + oidomain = idomain.reshape((nlay, nrow * ncol)) + for mkey in (0, 1): + ml = new_sim.get_model(f"ml_{mkey}") + grid = ml.modelgrid + nodes = mfsplit._grid_info[mkey][-1] + ncpl += grid.ncpl + + # a stack that is inactive or passes flow vertically in one layer + # keeps its idomain + for value in (0, -1): + nexpected = np.count_nonzero(oidomain[:, nodes] == value) + if np.count_nonzero(grid.idomain == value) != nexpected: + raise AssertionError( + f"Model {mkey} idomain {value} cells were not remapped" + ) + + heads[mkey] = ml.output.head().get_alldata()[-1] + + # a stack that is inactive in every layer is excluded + nactive = np.count_nonzero(np.any(oidomain != 0, axis=0)) + if ncpl != nactive: + raise AssertionError("Split models do not include all active cells") + + new_heads = mfsplit.reconstruct_array(heads) + + idx = idomain > 0 + err_msg = "Heads from original and split models do not match" + np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg) + + +def test_structured_to_disv_reconstruct_recarray(): + sim_path = get_example_data_path() / "mf6-freyberg" + + sim = MFSimulation.load(sim_ws=sim_path) + gwf = sim.get_model() + modelgrid = gwf.modelgrid + + array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int) + array[modelgrid.nrow // 2 :, :] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array, to_disv=True) + + # riv and wel are split across both models, chd is in one model + for pkgtype in ("riv", "wel", "chd"): + original = getattr(gwf, pkgtype).stress_period_data.get_data(0) + recarrays = {} + for mkey in (0, 1): + ml = new_sim.get_model(f"freyberg_{mkey}") + if hasattr(ml, pkgtype): + recarrays[mkey] = getattr(ml, pkgtype).stress_period_data.get_data(0) + + new_recarray = mfsplit.reconstruct_recarray(recarrays) + + # cellids are returned as the layer, row, column of the original model + onodes = modelgrid.get_node([tuple(cid) for cid in original.cellid]) + nnodes = modelgrid.get_node([tuple(cid) for cid in new_recarray.cellid]) + oidx = np.argsort(onodes) + nidx = np.argsort(nnodes) + + err_msg = f"Reconstructed {pkgtype} recarray does not match the original" + if not np.array_equal(np.array(onodes)[oidx], np.array(nnodes)[nidx]): + raise AssertionError(err_msg) + + for name in original.dtype.names: + if name == "cellid": + continue + np.testing.assert_allclose( + new_recarray[name][nidx], original[name][oidx], err_msg=err_msg + ) + + +def test_vertex_reconstruct_recarray(): + sim_path = get_example_data_path() / "mf6" / "test003_gwftri_disv" + + sim = MFSimulation.load(sim_ws=sim_path) + gwf = sim.get_model() + modelgrid = gwf.modelgrid + + array = np.zeros((modelgrid.ncpl,), dtype=int) + array[0:85] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array) + + for pname in ("chd_left", "chd_right"): + original = gwf.get_package(pname).stress_period_data.get_data(0) + recarrays = {} + for mkey in (0, 1): + ml = new_sim.get_model(f"gwf_1_{mkey}") + pkg = ml.get_package(pname) + if pkg is not None: + recarrays[mkey] = pkg.stress_period_data.get_data(0) + + new_recarray = mfsplit.reconstruct_recarray(recarrays) + + # cellids are returned as the layer, node of the original model + onodes = sorted(cid[-1] for cid in original.cellid) + nnodes = sorted(cid[-1] for cid in new_recarray.cellid) + + err_msg = f"Reconstructed {pname} recarray does not match the original" + if onodes != nnodes: + raise AssertionError(err_msg) + + +def test_structured_to_disv_model_geometry(): + sim_path = get_example_data_path() / "mf6-freyberg" + + sim = MFSimulation.load(sim_ws=sim_path) + gwf = sim.get_model() + gwf.modelgrid.set_coord_info(xoff=1000.0, yoff=2000.0, angrot=15.0) + modelgrid = gwf.modelgrid + + array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int) + array[:, modelgrid.ncol // 2 :] = 1 + + mfsplit = Mf6Splitter(sim) + mfsplit.split_model(array, to_disv=True) + + xcenters = modelgrid.xcellcenters.ravel() + ycenters = modelgrid.ycellcenters.ravel() + areas = modelgrid.area.ravel() + for mkey, model in mfsplit._model_dict.items(): + grid = model.modelgrid + nodes = mfsplit._grid_info[mkey][-1] + + err_msg = f"Model {mkey} cells are not coincident with original cells" + np.testing.assert_allclose(grid.xcellcenters, xcenters[nodes], err_msg=err_msg) + np.testing.assert_allclose(grid.ycellcenters, ycenters[nodes], err_msg=err_msg) + + # MODFLOW 6 calculates a CELL2D area that is positive only when the + # vertices are listed in clockwise order (Disv.f90 get_cell2d_area) + cell2d = model.disv.cell2d.array + vertices = model.disv.vertices.array + icvert = np.column_stack([cell2d[f"icvert_{iv}"] for iv in range(4)]) + x = vertices["xv"][icvert] + y = vertices["yv"][icvert] + cell_area = -0.5 * np.sum( + x * np.roll(y, -1, axis=1) - np.roll(x, -1, axis=1) * y, axis=1 + ) + err_msg = f"Model {mkey} CELL2D vertices are not in clockwise order" + np.testing.assert_allclose(cell_area, areas[nodes], err_msg=err_msg) + + +def test_to_disv_grid_types(): + sim_path = get_example_data_path() / "mf6" / "test003_gwftri_disv" + sim = MFSimulation.load(sim_ws=sim_path) + modelgrid = sim.get_model().modelgrid + + array = np.zeros((modelgrid.ncpl,), dtype=int) + array[0:85] = 1 + + mfsplit = Mf6Splitter(sim) + with pytest.warns(UserWarning, match="already a DISV model"): + mfsplit.split_model(array, to_disv=True) + + sim_path = get_example_data_path() / "mf6" / "test006_gwf3" + sim = MFSimulation.load(sim_ws=sim_path) + modelgrid = sim.get_model().modelgrid + + array = np.zeros((modelgrid.nnodes,), dtype=int) + array[65:] = 1 + + mfsplit = Mf6Splitter(sim) + with pytest.raises(ValueError, match="DISU"): + mfsplit.split_model(array, to_disv=True) + + +@requires_exe("mf6") +@requires_pkg("h5py") +def test_save_load_node_mapping_structured_to_disv(function_tmpdir): + sim_path = get_example_data_path() / "mf6-freyberg" + new_sim_path = function_tmpdir / "split_model" + hdf_file = new_sim_path / "node_map.hdf5" + + sim = MFSimulation.load(sim_ws=sim_path) + sim.set_sim_path(function_tmpdir) + sim.write_simulation() + sim.run_simulation() + + gwf = sim.get_model() + modelgrid = gwf.modelgrid + original_heads = gwf.output.head().get_alldata()[-1] + + array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int) + array[modelgrid.nrow // 2 :, :] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array, to_disv=True) + new_sim.set_sim_path(new_sim_path) + new_sim.write_simulation() + new_sim.run_simulation() + original_node_map = mfsplit._node_map + + mfsplit.save_node_mapping(hdf_file) + + new_sim2 = MFSimulation.load(sim_ws=new_sim_path) + mfsplit2 = Mf6Splitter.load_node_mapping(hdf_file) + + if original_node_map != mfsplit2._node_map: + raise AssertionError("Node map read/write not returning proper values") + + array_dict = {} + for mkey in (0, 1): + ml = new_sim2.get_model(f"freyberg_{mkey}") + if ml.modelgrid.grid_type != "vertex": + raise AssertionError("Loaded node mapping grid type is not vertex") + array_dict[mkey] = ml.output.head().get_alldata()[-1] + + new_heads = mfsplit2.reconstruct_array(array_dict) + + idx = modelgrid.idomain != 0 + err_msg = "Heads from original and split models do not match" + np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg) + + @requires_exe("mf6") @pytest.mark.slow def test_model_with_lak_sfr_mvr(function_tmpdir): @@ -174,6 +489,144 @@ def test_model_with_lak_sfr_mvr(function_tmpdir): np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg) +@requires_exe("mf6") +def test_structured_to_disv_multi_model(function_tmpdir): + nlay, nrow, ncol = 1, 10, 10 + idomain = np.ones((nlay, nrow, ncol), dtype=int) + idomain[:, 0:2, 0:2] = 0 + + sim = flopy.mf6.MFSimulation(sim_name="mm", sim_ws=function_tmpdir, exe_name="mf6") + flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(100.0, 10, 1.0)]) + ims_gwf = flopy.mf6.ModflowIms(sim, complexity="simple", filename="gwf.ims") + ims_gwt = flopy.mf6.ModflowIms( + sim, + complexity="simple", + linear_acceleration="bicgstab", + filename="gwt.ims", + ) + + dis_kwargs = { + "nlay": nlay, + "nrow": nrow, + "ncol": ncol, + "delr": 100.0, + "delc": 100.0, + "top": 10.0, + "botm": [0.0], + "idomain": idomain, + } + + gwf = flopy.mf6.ModflowGwf(sim, modelname="gwf", save_flows=True) + flopy.mf6.ModflowGwfdis(gwf, **dis_kwargs) + flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True, k=1.0) + flopy.mf6.ModflowGwfic(gwf, strt=10.0) + flopy.mf6.ModflowGwfsto(gwf, ss=1e-5, iconvert=0) + chd = [[(0, i, 2), 10.0, 1.0] for i in range(2, nrow)] + chd += [[(0, i, ncol - 1), 9.0, 0.0] for i in range(nrow)] + flopy.mf6.ModflowGwfchd( + gwf, pname="chd-1", auxiliary=["conc"], stress_period_data=chd + ) + flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord="gwf.hds", + budget_filerecord="gwf.cbc", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + gwt = flopy.mf6.ModflowGwt(sim, modelname="gwt") + flopy.mf6.ModflowGwtdis(gwt, **dis_kwargs) + flopy.mf6.ModflowGwtic(gwt, strt=0.0) + flopy.mf6.ModflowGwtmst(gwt, porosity=0.2) + flopy.mf6.ModflowGwtadv(gwt, scheme="upstream") + flopy.mf6.ModflowGwtssm(gwt, sources=[["chd-1", "AUX", "conc"]]) + flopy.mf6.ModflowGwtoc( + gwt, + concentration_filerecord="gwt.ucn", + saverecord=[("CONCENTRATION", "ALL")], + ) + flopy.mf6.ModflowGwfgwt(sim, exgmnamea="gwf", exgmnameb="gwt") + sim.register_ims_package(ims_gwf, ["gwf"]) + sim.register_ims_package(ims_gwt, ["gwt"]) + sim.write_simulation() + sim.run_simulation() + + original_conc = gwt.output.concentration().get_alldata()[-1] + + array = np.zeros((nrow, ncol), dtype=int) + array[:, ncol // 2 :] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_multi_model(array, to_disv=True) + new_sim.set_sim_path(function_tmpdir / "split_model") + new_sim.write_simulation() + new_sim.run_simulation() + + conc = {} + for mkey in (0, 1): + for mname in (f"gwf_{mkey}", f"gwt_{mkey}"): + if new_sim.get_model(mname).modelgrid.grid_type != "vertex": + raise AssertionError(f"Model {mname} is not a DISV model") + + conc[mkey] = new_sim.get_model(f"gwt_{mkey}").output.concentration() + conc[mkey] = conc[mkey].get_alldata()[-1] + + new_conc = mfsplit.reconstruct_array(conc) + + idx = idomain != 0 + err_msg = "Concentrations from original and split models do not match" + np.testing.assert_allclose( + new_conc[idx], original_conc[idx], atol=1e-6, err_msg=err_msg + ) + + +@requires_exe("mf6") +@pytest.mark.slow +def test_structured_to_disv_with_lak_sfr_mvr(function_tmpdir): + sim_path = get_example_data_path() / "mf6" / "test045_lake2tr" + + sim = MFSimulation.load(sim_ws=sim_path) + sim.set_sim_path(function_tmpdir) + sim.write_simulation() + sim.run_simulation() + + gwf = sim.get_model() + modelgrid = gwf.modelgrid + + array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int) + array[0:14, :] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array, to_disv=True) + + new_sim.set_sim_path(function_tmpdir / "split_model") + new_sim.write_simulation() + new_sim.run_simulation() + + original_heads = gwf.output.head().get_alldata()[-1] + + ml0 = new_sim.get_model("lakeex2a_0") + ml1 = new_sim.get_model("lakeex2a_1") + for ml in (ml0, ml1): + if ml.modelgrid.grid_type != "vertex": + raise AssertionError("Split model is not a DISV model") + + # the advanced packages remap cellids through the new grid type + for pkgtype in ("lak", "sfr", "mvr", "evta", "rcha"): + if ml.get_package(pkgtype) is None: + raise AssertionError(f"{pkgtype} package was not split") + + heads0 = ml0.output.head().get_alldata()[-1] + heads1 = ml1.output.head().get_alldata()[-1] + + new_heads = mfsplit.reconstruct_array({0: heads0, 1: heads1}) + + idx = modelgrid.idomain != 0 + err_msg = "Heads from original and split models do not match" + np.testing.assert_allclose( + new_heads[idx], original_heads[idx], atol=1e-4, err_msg=err_msg + ) + + @requires_exe("mf6") @requires_pkg("pymetis") @pytest.mark.slow @@ -659,7 +1112,8 @@ def test_empty_ssm(function_tmpdir): @requires_exe("mf6") -def test_transient_array(function_tmpdir): +@pytest.mark.parametrize("to_disv", [False, True]) +def test_transient_array(function_tmpdir, to_disv): name = "tarr" new_sim_path = function_tmpdir / f"{name}_split_model" nper = 3 @@ -737,7 +1191,7 @@ def test_transient_array(function_tmpdir): sarr = np.ones((nrow, ncol), dtype=int) sarr[:, int(ncol / 2) :] = 2 mfsplit = Mf6Splitter(sim) - new_sim = mfsplit.split_model(sarr) + new_sim = mfsplit.split_model(sarr, to_disv=to_disv) for name in new_sim.model_names: g = new_sim.get_model(name) diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 690cf306e..077017073 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -2001,21 +2001,43 @@ def get_plottable_layer_array(self, a, layer): assert plotarray.shape == required_shape, msg return plotarray + def get_cell_iverts(self, nodes=None): + """ + Get the vertex numbers that define one or more model cells + + Parameters + ---------- + nodes : int, list, np.ndarray + optional two-dimensional node numbers. Vertex numbers are returned + for every cell in a layer when nodes is None. + + Returns + ------- + np.ndarray : array of vertex numbers with shape (nnodes, 4). The + vertices of a cell are ordered clockwise from its upper left + corner. + + """ + if nodes is None: + nodes = np.arange(self.ncpl, dtype=int) + else: + nodes = np.ravel(np.asarray(nodes, dtype=int)) + + i, j = np.divmod(nodes, self.ncol) + iverts = np.empty((4, nodes.size), dtype=int) + iverts[0] = i * (self.ncol + 1) + j + iverts[1] = iverts[0] + 1 + iverts[2] = iverts[1] + self.ncol + 1 + iverts[3] = iverts[0] + self.ncol + 1 + return iverts.T + def _set_structured_iverts(self): """ Build a list of the vertices that define each model cell and the x, y pair for each vertex """ - rowarr = np.repeat(np.arange(self.nrow, dtype=int), self.ncol) - colarr = np.tile(np.arange(self.ncol, dtype=int), self.nrow) - - iverts = np.empty((4, self.ncpl), dtype=int) - iverts[0] = rowarr * (self.ncol + 1) + colarr - iverts[1] = rowarr * (self.ncol + 1) + colarr + 1 - iverts[2] = (rowarr + 1) * (self.ncol + 1) + colarr + 1 - iverts[3] = (rowarr + 1) * (self.ncol + 1) + colarr - self._iverts = iverts.T.tolist() + self._iverts = self.get_cell_iverts().tolist() return def _build_structured_iverts(self, i, j): diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index 8b333cb5e..137b041e2 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -1,4 +1,5 @@ import inspect +import warnings import numpy as np @@ -153,6 +154,8 @@ def __init__(self, sim, modelname=None): self._ncpl = self._modelgrid.nnodes self._shape = self._modelgrid.shape self._grid_type = self._modelgrid.grid_type + self._new_grid_type = self._grid_type + self._to_disv = False self._node_map = {} self._node_map_r = {} self._new_connections = None @@ -283,10 +286,12 @@ def _map_grids_to_hdf5(self, f): } mgg = f.create_group("modelgrids") - grid_type = f["grid_type"][0].decode("utf8") mkeys = f["mkey"][:] mnames = [self._modelname] + list(f["new_modelnames"][:]) grids = [self._modelgrid, ] + [self._model_dict[mk].modelgrid for mk in mkeys] + grid_types = [f["grid_type"][0].decode("utf8")] + [ + f["new_grid_type"][0].decode("utf8") + ] * len(mkeys) for ix, name in enumerate(mnames): if hasattr(name, "decode"): @@ -295,6 +300,7 @@ def _map_grids_to_hdf5(self, f): # name = name.decode("utf8") grd_grp = mgg.create_group(name) grid = grids[ix] + grid_type = grid_types[ix] if grid_type == "structured": delc = grid.delc dc = grd_grp.create_dataset("delc", (len(delc),), dtype=float, **kwargs) @@ -411,6 +417,9 @@ def save_node_mapping(self, filename): gt_ds = f.create_dataset("grid_type", (1,), dtype=string_dt) gt_ds[:] = [self._grid_type,] + ngt_ds = f.create_dataset("new_grid_type", (1,), dtype=string_dt) + ngt_ds[:] = [self._new_grid_type,] + mname_ds = f.create_dataset("modelname", (1,), dtype=string_dt) mname_ds[:] = [self._modelname,] @@ -561,10 +570,15 @@ def construct_modelgrid(f, name, grid_type): # construct representation of original model geometry modelname = f["modelname"][0].decode("utf8") grid_type = f["grid_type"][0].decode("utf8") + new_grid_type = grid_type + if "new_grid_type" in f: + new_grid_type = f["new_grid_type"][0].decode("utf8") + grid = construct_modelgrid(f, modelname, grid_type) ml = FakeModel(modelname, grid) sim = FakeSim({modelname: ml}) mfs = Mf6Splitter(sim) + mfs._new_grid_type = new_grid_type mfs._model_dict = {} mfs._new_ncpl = {} @@ -573,10 +587,10 @@ def construct_modelgrid(f, name, grid_type): modelnames = [i.decode("utf8") for i in f["new_modelnames"][:]] for ix, mname in enumerate(modelnames): mkey = mkeys[ix] - sgrid = construct_modelgrid(f, mname, grid_type) + sgrid = construct_modelgrid(f, mname, new_grid_type) sml = FakeModel(mname, sgrid) mfs._model_dict[mkey] = sml - if grid_type in ("structured", "vertex"): + if new_grid_type in ("structured", "vertex"): mfs._new_ncpl[mkey] = sgrid.ncpl else: mfs._new_ncpl[mkey] = sgrid.nnodes @@ -584,8 +598,9 @@ def construct_modelgrid(f, name, grid_type): f.close() # create additional splitting data efficient reconstruction - split_array = np.zeros((mfs._ncpl,), dtype=int) - model_array = np.zeros((mfs._ncpl,), dtype=int) + # cells that were not remapped keep a model number of -1 + split_array = np.full((mfs._ncpl,), -1, dtype=int) + model_array = np.full((mfs._ncpl,), -1, dtype=int) for k, v in node_map.items(): k = int(k) model_array[k] = v[0] @@ -867,9 +882,9 @@ def reconstruct_recarray(self, recarrays): remapper = self.reversed_node_map[mkey] orec = recarray.copy() modelgrid = self._model_dict[mkey].modelgrid - if self._grid_type in ("structured", "vertex"): + if self._new_grid_type in ("structured", "vertex"): layer = [i[0] for i in orec.cellid] - if self._grid_type == "structured": + if self._new_grid_type == "structured": cellid = [(0, i[1], i[2]) for i in orec.cellid] node = modelgrid.get_node(cellid) else: @@ -879,7 +894,7 @@ def reconstruct_recarray(self, recarrays): new_node = [remapper[i] for i in node if i in remapper] - if modelgrid.grid_type == "structured": + if self._modelgrid.grid_type == "structured": if self._modelgrid is None: new_cellid = list( zip(*np.unravel_index(new_node, self._shape)) @@ -889,7 +904,7 @@ def reconstruct_recarray(self, recarrays): new_cellid = [ (layer[ix], i[1], i[2]) for ix, i in enumerate(new_cellid) ] - elif modelgrid.grid_type == "vertex": + elif self._modelgrid.grid_type == "vertex": new_cellid = [(layer[ix], i) for ix, i in enumerate(new_node)] else: new_cellid = [(i,) for i in new_node] @@ -1013,7 +1028,7 @@ def _remap_nodes(self, array): ) grid_info = {} - if self._modelgrid.grid_type == "structured": + if self._modelgrid.grid_type == "structured" and not self._to_disv: a = array.reshape(self._modelgrid.nrow, self._modelgrid.ncol) for m in np.unique(a): cells = np.asarray(a == m).nonzero() @@ -1038,8 +1053,16 @@ def _remap_nodes(self, array): np.ravel(mapping), ] else: + # a DISV model has no use for a stack of cells that is inactive + # in every layer, so those cells are not carried into it + mask = np.ones((self._ncpl,), dtype=bool) + if self._to_disv: + for arr in idomain: + mask &= arr == 0 + mask = ~mask + for m in mkeys: - cells = np.asarray(array == m).nonzero()[0] + cells = np.asarray((array == m) & mask).nonzero()[0] mapping = np.zeros((len(cells),), dtype=int) mapping[:] = cells grid_info[m] = [(len(cells),), None, None, mapping] @@ -1075,8 +1098,14 @@ def _remap_nodes(self, array): exchange_meta = {i: {} for i in mkeys} usg_meta = {i: {} for i in mkeys} for node, conn in self._connection.items(): + if node not in self._node_map: + continue + mdl, nnode = self._node_map[node] for ix, cnode in enumerate(conn): + if cnode not in self._node_map: + continue + cmdl, cnnode = self._node_map[cnode] if cmdl == mdl: if nnode in new_connections[mdl]["internal"]: @@ -1327,6 +1356,58 @@ def _remap_cell2d(self, item, cell2d, mapped_data): return mapped_data + def _structured_to_disv(self, mapped_data): + """ + Method to create the DISV grid of models split from a structured model + + Parameters + ---------- + mapped_data : dict + dictionary of remapped package data + + Returns + ------- + dict + """ + modelgrid = self._modelgrid + ncol = modelgrid.ncol + xedge, yedge = modelgrid.xyedges + xcenter, ycenter = modelgrid.xycenters + vert_dtype = np.dtype([("iv", int), ("xv", float), ("yv", float)]) + cell_dtype = np.dtype( + [("icell2d", int), ("xc", float), ("yc", float), ("ncvert", int)] + + [(f"icvert_{iv}", int) for iv in range(4)] + ) + + for mkey in self._model_dict.keys(): + nodes = self._grid_info[mkey][-1] + # corners shared by cells in the model collapse to one vertex + iverts = modelgrid.get_cell_iverts(nodes) + overts, iverts = np.unique(iverts, return_inverse=True) + iverts = iverts.reshape((nodes.size, 4)) + + i, j = np.divmod(overts, ncol + 1) + vertices = np.recarray((overts.size,), dtype=vert_dtype) + vertices["iv"] = np.arange(overts.size, dtype=int) + vertices["xv"] = xedge[j] + vertices["yv"] = yedge[i] + + i, j = np.divmod(nodes, ncol) + cell2d = np.recarray((nodes.size,), dtype=cell_dtype) + cell2d["icell2d"] = np.arange(nodes.size, dtype=int) + cell2d["xc"] = xcenter[j] + cell2d["yc"] = ycenter[i] + cell2d["ncvert"] = 4 + for iv in range(4): + cell2d[f"icvert_{iv}"] = iverts[:, iv] + + mapped_data[mkey]["ncpl"] = self._new_ncpl[mkey] + mapped_data[mkey]["nvert"] = overts.size + mapped_data[mkey]["vertices"] = vertices + mapped_data[mkey]["cell2d"] = cell2d + + return mapped_data + def _remap_filerecords(self, item, value, mapped_data, namfile=False): """ Method to create new file record names and map them to their @@ -3133,10 +3214,10 @@ def _new_node_to_cellid(self, model, new_node, layers, idx): """ new_node = new_node[idx].astype(int) - if self._modelgrid.grid_type == "structured": + if self._new_grid_type == "structured": new_node += layers[idx] * model.modelgrid.ncpl new_cellids = model.modelgrid.get_lrc(new_node.astype(int)) - elif self._modelgrid.grid_type == "vertex": + elif self._new_grid_type == "vertex": new_cellids = [tuple(cid) for cid in zip(layers[idx], new_node)] else: @@ -3279,6 +3360,9 @@ def _remap_package(self, package, ismvr=False): continue if item in ("delr", "delc"): + if self._to_disv: + continue + for mkey, d in self._grid_info.items(): if item == "delr": i0, i1 = d[2] @@ -3288,6 +3372,9 @@ def _remap_package(self, package, ismvr=False): mapped_data[mkey][item] = value.array[i0 : i1 + 1] elif item in ("nrow", "ncol"): + if self._to_disv: + continue + for mkey, d in self._grid_info.items(): if item == "nrow": i0, i1 = d[1] @@ -3330,6 +3417,9 @@ def _remap_package(self, package, ismvr=False): elif isinstance(value, mfdataarray.MFArray): mapped_data = self._remap_array(item, value, mapped_data) + if self._to_disv: + mapped_data = self._structured_to_disv(mapped_data) + elif isinstance(package, modflow.ModflowGwfhfb): mapped_data = self._remap_hfb(package, mapped_data) @@ -3550,8 +3640,14 @@ def _remap_package(self, package, ismvr=False): elif isinstance(value, mfdatalist.MFList): mapped_data[mkey][item] = value.array + package_type = package.package_type + pname = package.name[0] + if self._to_disv and package_type == "dis": + package_type = "disv" + pname = None + pak_cls = PackageContainer.package_factory( - package.package_type, self._model_type + package_type, self._model_type ) paks = {} for mdl, data in mapped_data.items(): @@ -3570,7 +3666,7 @@ def _remap_package(self, package, ismvr=False): self._new_sim.simulation_data.max_columns_of_data = max_cols paks[mdl] = pak_cls( - self._model_dict[mdl], pname=package.name[0], **data + self._model_dict[mdl], pname=pname, **data ) if observations is not None: @@ -3723,7 +3819,7 @@ def _create_exchanges(self): node1 = exg[1] for layer in range(self._modelgrid.nlay): - if self._modelgrid.grid_type == "structured": + if self._new_grid_type == "structured": tmpnode0 = node0 + (ncpl0 * layer) tmpnode1 = node1 + (ncpl1 * layer) cellidm0 = modelgrid0.get_lrc([tmpnode0])[ @@ -3732,7 +3828,7 @@ def _create_exchanges(self): cellidm1 = modelgrid1.get_lrc([tmpnode1])[ 0 ] - elif self._modelgrid.grid_type == "vertex": + elif self._new_grid_type == "vertex": cellidm0 = (layer, node0) cellidm1 = (layer, node1) else: @@ -3892,7 +3988,37 @@ def create_multi_model_exchanges(self, mname0, mname1): filename=filename, ) - def split_model(self, array, sim_ws=None): + def _set_new_grid_type(self, to_disv): + """ + Method to set the grid type of the models created by splitting + + Parameters + ---------- + to_disv : bool + flag to write DISV models when splitting a structured model + + """ + grid_type = self._modelgrid.grid_type + self._to_disv = False + self._new_grid_type = grid_type + if not to_disv: + return + + if grid_type == "unstructured": + raise ValueError("DISV models cannot be split from a DISU model") + + if grid_type == "vertex": + warnings.warn( + "the model is already a DISV model, to_disv is ignored", + UserWarning, + stacklevel=3, + ) + return + + self._to_disv = True + self._new_grid_type = "vertex" + + def split_model(self, array, sim_ws=None, to_disv=False): """ User method to split a model based on an array @@ -3906,6 +4032,11 @@ def split_model(self, array, sim_ws=None): optional directory path for writing the new simulation to. This parameter is recommended when the model contains external files and the user would like to preserve external linkages while splitting. + to_disv : bool + optional flag to write DISV models when splitting a structured model. + A DISV model carries only the cells assigned to it, not the inactive + cells that pad a structured model to its bounding box. Stacks of cells + that are inactive in every layer are excluded. Default is False. Returns ------- @@ -3917,6 +4048,8 @@ def split_model(self, array, sim_ws=None): "is part of a split simulation" ) + self._set_new_grid_type(to_disv) + if sim_ws is None: self._keep_external = False sim_ws = self._sim.sim_path @@ -3970,7 +4103,7 @@ def split_model(self, array, sim_ws=None): return self._new_sim - def split_multi_model(self, array, sim_ws=None): + def split_multi_model(self, array, sim_ws=None, to_disv=False): """ Method to split integrated models such as GWF-GWT or GWF-GWE models. Note: this method will not work to split multiple connected GWF models @@ -3985,6 +4118,11 @@ def split_multi_model(self, array, sim_ws=None): optional directory path for writing the new simulation to. This parameter is recommended when the model contains external files and the user would like to preserve external linkages while splitting. + to_disv : bool + optional flag to write DISV models when splitting a structured model. + A DISV model carries only the cells assigned to it, not the inactive + cells that pad a structured model to its bounding box. Stacks of cells + that are inactive in every layer are excluded. Default is False. Returns ------- @@ -4048,10 +4186,10 @@ def split_multi_model(self, array, sim_ws=None): int(i): f"{gwf_base}_{i}" for i in model_labels } - new_sim = self.split_model(array) + new_sim = self.split_model(array, to_disv=to_disv) for mname in model_names[1:]: self.switch_models(modelname=mname, remap_nodes=False) - new_sim = self.split_model(array, sim_ws=sim_ws) + new_sim = self.split_model(array, sim_ws=sim_ws, to_disv=to_disv) for mbase in model_names[1:]: for label in model_labels: