diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 4246d0513..10df3ff02 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1981,3 +1981,134 @@ def test_reconstruct_recarray(function_tmpdir): vrec = vrecarray[ix] for name in recarray.dtype.names: assert rec[name] == vrec[name], "Recarray reconstruction failed" + + +@requires_exe("mf6") +def test_fjf_obs_structured(function_tmpdir): + sim_ws = function_tmpdir + split_ws = function_tmpdir / "split_model" + + nlay = 1 + nrow = 10 + ncol = 10 + delc = np.full((nrow,), 10) + delr = np.full((ncol,), 10) + top = np.full((nrow, ncol), 50) + botm = np.full((nlay, nrow, ncol), 0) + idomain = np.ones(botm.shape, dtype=int) + + sim = flopy.mf6.MFSimulation(sim_ws=sim_ws) + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + gwf = flopy.mf6.ModflowGwf(sim) + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=idomain, + ) + npf = flopy.mf6.ModflowGwfnpf(gwf, k=5, k33=1) + ic = flopy.mf6.ModflowGwfic(gwf, strt=top - 5) + + chd_rec = [(0, i, 0, 45) for i in range(nrow)] + ghb_rec = [(0, i, 9, 35, 25) for i in range(nrow)] + chd = flopy.mf6.ModflowGwfchd( + gwf, stress_period_data={0: chd_rec}, pname="chd_left" + ) + ghb = flopy.mf6.ModflowGwfghb( + gwf, stress_period_data={0: ghb_rec}, pname="ghb_right" + ) + + obs = flopy.mf6.ModflowUtlobs( + gwf, + continuous={ + "obs.csv": [ + ("faceflow1", "flow-ja-face", (0, 3, 3), (0, 3, 4)), + ("h1", "head", (0, 3, 4)), + ("faceflow2", "flow-ja-face", (0, 6, 6), (0, 6, 7)), + ("h2", "head", (0, 6, 7)), + ] + }, + ) + + sim.write_simulation() + sim.run_simulation() + + array = np.zeros(top.shape, dtype=int) + array[:, (ncol // 2) :] = 1 + mfs = flopy.mf6.utils.Mf6Splitter(sim) + new_sim = mfs.split_model(array, sim_ws=split_ws) + new_sim.write_simulation() + success, _ = new_sim.run_simulation() + assert success, "split model run failed, obs remapping may have failed" + + for mdl in range(2): + obs_ra = new_sim.get_model(f"{gwf.name}_{mdl}").obs.continuous.data[ + f"obs_{mdl}.csv" + ] + assert len(obs_ra) == 2, "number remapped of observation records incorrect" + id2 = obs_ra[obs_ra.obsname == f"faceflow{mdl + 1}"].id2[0] + assert isinstance(id2, tuple), "obs id2 datatype is incorrect" + + if mdl == 0: + v_ra = obs.continuous.data["obs.csv"] + v_id = v_ra[v_ra.obsname == f"faceflow{mdl + 1}"].id2[0] + + assert id2 == v_id, "id2 cellid not properly set" + + +@requires_exe("mf6") +def test_fjf_obs_vertex(function_tmpdir): + sim_ws = get_example_data_path() / "mf6" / "test003_gwfs_disv" + split_ws = function_tmpdir / "split_model" + + sim = flopy.mf6.MFSimulation.load(sim_ws=sim_ws) + sim.set_sim_path(function_tmpdir) + + gwf = sim.get_model() + modelgrid = gwf.modelgrid + ncpl = modelgrid.ncpl + + obs = flopy.mf6.ModflowUtlobs( + gwf, + continuous={ + "obs.csv": [ + ("faceflow1", "flow-ja-face", (0, 24), (0, 34)), + ("h1", "head", (0, 34)), + ("faceflow2", "flow-ja-face", (0, 64), (0, 74)), + ("h2", "head", (0, 74)), + ] + }, + ) + + sim.write_simulation() + sim.run_simulation() + + array = np.zeros((ncpl,), dtype=int) + array[(ncpl // 2) :] = 1 + + mfs = flopy.mf6.utils.Mf6Splitter(sim) + new_sim = mfs.split_model(array, sim_ws=split_ws) + new_sim.write_simulation() + success, _ = new_sim.run_simulation() + assert success, "split model run failed, obs remapping may have failed" + + for mdl in range(2): + obs_ra = new_sim.get_model(f"{gwf.name}_{mdl}").obs.continuous.data[ + f"obs_{mdl}.csv" + ] + assert len(obs_ra) == 2, "number remapped of observation records incorrect" + id2 = obs_ra[obs_ra.obsname == f"faceflow{mdl + 1}"].id2[0] + assert isinstance(id2, tuple), "obs id2 datatype is incorrect" + + if mdl == 0: + v_ra = obs.continuous.data["obs.csv"] + v_id = v_ra[v_ra.obsname == f"faceflow{mdl + 1}"].id2[0] + + assert id2 == v_id, "id2 cellid not properly set" diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index d2676a712..07ddd7512 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -2929,10 +2929,10 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): "structured", "vertex", ): - layers2 = [ + layers2 = np.array([ cid[0] if cid is not None else None for cid in cellid2 - ] + ], dtype=object) if self._modelgrid.grid_type == "structured": cellid2 = [ ( @@ -2949,7 +2949,7 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): ] node2 = self._modelgrid.get_node( - list(cellid2[conv_idx]) + [cellid2[cv_ix] for cv_ix in conv_idx] ) new_node2 = np.full( (len(recarray),), None, dtype=object @@ -2964,10 +2964,10 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): if ix in conv_idx: continue else: - new_node2.append(None) - new_model2.append(new_model1[ix]) + new_node2[ix] = None + new_model2[ix] = int(new_model1[ix]) - if not np.allclose(new_model1, new_model2): + if not np.allclose(new_model1, new_model2.astype(int)): raise AssertionError( "One or more observation records cross model boundaries" ) @@ -2978,20 +2978,11 @@ def _remap_obs(self, package, mapped_data, remapper, pkg_type=None): for mkey, model in self._model_dict.items(): idx = np.asarray(new_model2 == mkey).nonzero() tmp_node = new_node2[idx] - cidx = np.asarray(tmp_node != None).nonzero() # noqa: E711 - tmp_cellid = model.modelgrid.get_lrc( - tmp_node[cidx].to_list() + tmp_layers = layers2[idx] + cidx = np.asarray(tmp_node != None).nonzero() # noqa: E711 + tmp_cellid = self._new_node_to_cellid( + model, tmp_node, tmp_layers, cidx ) - if self._modelgrid.grid_type in ( - "structured", - "vertex", - ): - tmp_layers = layers2[cidx] - tmp_cellid = [ - (tmp_layers[ix],) + cid[1:] - for ix, cid in enumerate(tmp_cellid) - ] - tmp_node[cidx] = tmp_cellid new_cellid2[idx] = tmp_node else: @@ -3169,7 +3160,7 @@ def _new_node_to_cellid(self, model, new_node, layers, idx): new_node = new_node[idx].astype(int) if self._modelgrid.grid_type == "structured": - new_node += layers[idx] * model.modelgrid.ncpl + new_node += layers[idx].astype(int) * model.modelgrid.ncpl new_cellids = model.modelgrid.get_lrc(new_node.astype(int)) elif self._modelgrid.grid_type == "vertex": new_cellids = [tuple(cid) for cid in zip(layers[idx], new_node)]