Fix(_remap_obs): fix support for GWF observation types with 2 cellids

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"

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)]