fix(mf6): fix external binary files for vertex grids (#1839)

autotest/test_mf6.py

@@ -81,6 +81,8 @@
     ZoneBudget6,
 )
 from flopy.utils.observationfile import CsvFile
+from flopy.utils.triangle import Triangle
+from flopy.utils.voronoi import VoronoiGrid
 
 pytestmark = pytest.mark.mf6
 
@@ -502,6 +504,191 @@ def test_subdir(function_tmpdir):
     ), "Something wrong with model external paths"
 
 
+@requires_exe("mf6")
+def test_binary_write(function_tmpdir):
+    nlay, nrow, ncol = 2, 1, 10
+    shape2d = (nrow, ncol)
+    shape3d = (nlay, nrow, ncol)
+
+    # create binary data structured
+    top_data = {
+        "filename": "top.bin",
+        "binary": True,
+        "iprn": 0,
+        "data": 10.0,
+    }
+    botm_data = {
+        "filename": "botm.bin",
+        "binary": True,
+        "iprn": 0,
+        "data": np.array([
+            np.full(shape2d, 4.0, dtype=float),
+            np.full(shape2d, 0.0, dtype=float),
+        ]),
+    }
+    strt_data = {
+        "filename": "strt.bin",
+        "binary": True,
+        "iprn": 0,
+        "data": np.full(shape3d, 10., dtype=float),
+    }
+    rch_data = {
+        0: {
+            "filename": "recharge.bin",
+            "binary": True,
+            "iprn": 0,
+            "data": 0.000001,
+        },
+    }
+    chd_data = [
+        (1, 0, 0, 10.0, 1.0, 100.0),
+        (1, 0, ncol - 1, 5.0, 0.0, 100.0),
+    ]
+    chd_data = {
+        0: {
+            "filename": "chd.bin",
+            "binary": True,
+            "iprn": 0,
+            "data": chd_data,
+        },
+    }
+
+    sim = MFSimulation(sim_ws=str(function_tmpdir))
+    ModflowTdis(sim)
+    ModflowIms(sim, complexity="simple")
+    gwf = ModflowGwf(sim, print_input=True)
+    ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=1.0,
+        delc=1.0,
+        top=top_data,
+        botm=botm_data,
+    )
+    ModflowGwfnpf(
+        gwf,
+        icelltype=1,
+    )
+    ModflowGwfic(
+        gwf,
+        strt=10.0,
+    )
+    ModflowGwfchd(
+        gwf,
+        auxiliary=["conc", "something"],
+        stress_period_data=chd_data,
+    )
+    ModflowGwfrcha(gwf, recharge=rch_data)
+
+    sim.write_simulation()
+    success, buff = sim.run_simulation()
+    assert success
+
+
+@requires_exe("mf6")
+def test_vor_binary_write(function_tmpdir):
+
+    # build voronoi grid
+    boundary = [(0.0, 0.0), (0.0, 1.0), (10.0, 1.0), (10.0, 0.0)]
+    triangle_ws = function_tmpdir / "triangle"
+    triangle_ws.mkdir(parents=True, exist_ok=True)
+
+    tri = Triangle(
+        angle=30,
+        maximum_area=1.0,
+        model_ws=triangle_ws,
+    )
+    tri.add_polygon(boundary)
+    tri.build(verbose=False)
+    vor = VoronoiGrid(tri)
+
+    # problem dimensions
+    nlay = 2
+    shape3d = (nlay, vor.ncpl)
+
+    # build binary data
+    top_data = {
+        "filename": "top.bin",
+        "binary": True,
+        "iprn": 0,
+        "data": 10.0,
+    }
+    botm_data = {
+        "filename": "botm.bin",
+        "binary": True,
+        "iprn": 0,
+        "data": np.array([
+            np.full(vor.ncpl, 4.0, dtype=float),
+            np.full(vor.ncpl, 0.0, dtype=float),
+        ]),
+    }
+    strt_data = {
+        "filename": "strt.bin",
+        "binary": True,
+        "iprn": 0,
+        "data": np.full(shape3d, 10., dtype=float),
+    }
+    rch_data = {
+        0: {
+            "filename": "recharge.bin",
+            "binary": True,
+            "iprn": 0,
+            "data": np.full(vor.ncpl, 0.000001, dtype=float), #0.000001,
+        },
+    }
+    chd_data = [
+        (1, 0, 10.0, 1.0, 100.0),
+        (1, 1, 10.0, 1.0, 100.0),
+        (1, 2, 5.0, 0.0, 100.0),
+        (1, 3, 5.0, 0.0, 100.0),
+    ]
+    chd_data = {
+        0: {
+            "filename": "chd.bin",
+            "binary": True,
+            "data": chd_data,
+        },
+    }
+
+    # build model
+    sim = MFSimulation(sim_ws=str(function_tmpdir))
+    ModflowTdis(sim)
+    ModflowIms(sim, complexity="simple")
+    gwf = ModflowGwf(sim, print_input=True)
+    flopy.mf6.ModflowGwfdisv(
+        gwf,
+        nlay=nlay,
+        ncpl=vor.ncpl,
+        nvert=vor.nverts,
+        vertices=vor.get_disv_gridprops()["vertices"],
+        cell2d=vor.get_disv_gridprops()["cell2d"],
+        top=top_data,
+        botm=botm_data,
+        xorigin=0.0,
+        yorigin=0.0,
+    )
+    ModflowGwfnpf(
+        gwf,
+        icelltype=1,
+    )
+    ModflowGwfic(
+        gwf,
+        strt=strt_data,
+    )
+    ModflowGwfrcha(gwf, recharge=rch_data)
+    ModflowGwfchd(
+        gwf,
+        auxiliary=["conc", "something"],
+        stress_period_data=chd_data,
+    )
+    sim.write_simulation()
+    success, buff = sim.run_simulation()
+    assert success
+
+
+
 def test_binary_read(function_tmpdir):
     test_ex_name = "binary_read"
     nlay = 3

autotest/test_model_splitter.py

@@ -120,6 +120,7 @@ def test_unstructured_model_splitter(function_tmpdir):
 
 
 @requires_exe("mf6")
+@pytest.mark.slow
 def test_model_with_lak_sfr_mvr(function_tmpdir):
     sim_path = get_example_data_path() / "mf6" / "test045_lake2tr"
 
@@ -155,6 +156,7 @@ def test_model_with_lak_sfr_mvr(function_tmpdir):
 
 
 @requires_exe("mf6")
+@pytest.mark.slow
 def test_metis_splitting_with_lak_sfr(function_tmpdir):
     sim_path = get_example_data_path() / "mf6" / "test045_lake2tr"
 
@@ -199,7 +201,7 @@ def test_metis_splitting_with_lak_sfr(function_tmpdir):
 @requires_exe("mf6")
 def test_save_load_node_mapping(function_tmpdir):
     sim_path = get_example_data_path() / "mf6-freyberg"
-    new_sim_path = sim_path / "split_model"
+    new_sim_path = function_tmpdir / "mf6-freyberg/split_model"
     json_file = new_sim_path / "node_map.json"
     nparts = 5
 

flopy/mf6/data/mffileaccess.py

@@ -294,8 +294,8 @@ def _get_header(
                     precision=precision,
                     text=text,
                     ncpl=modelgrid.ncpl,
+                    m2=1,
                     ilay=ilay,
-                    m3=1,
                     pertim=pertim,
                     totim=totim,
                     kstp=1,
@@ -1097,7 +1097,7 @@ def _get_header(self, modelgrid, precision):
     def _get_cell_header(self, modelgrid):
         if modelgrid.grid_type == "structured":
             return [("layer", np.int32), ("row", np.int32), ("col", np.int32)]
-        elif modelgrid.grid_type == "vertex_layered":
+        elif modelgrid.grid_type == "vertex":
             return [("layer", np.int32), ("ncpl", np.int32)]
         else:
             return [("nodes", np.int32)]

flopy/utils/datafile.py

@@ -95,8 +95,8 @@ def __init__(self, filetype=None, precision="single"):
                         ("totim", floattype),
                         ("text", "a16"),
                         ("ncpl", "i4"),
+                        ("m2", "i4"),
                         ("ilay", "i4"),
-                        ("m3", "i4"),
                     ]
                 )
             elif self.header_type == "vardisu":