From 8113a8856a61fdbb7a743d4687ddb3fb68d99e62 Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Mon, 26 Jun 2023 11:55:41 -0500 Subject: [PATCH] fix(mf6): fix external binary files for vertex grids (#1839) --- autotest/test_mf6.py | 187 ++++++++++++++++++++++++++++++++ autotest/test_model_splitter.py | 4 +- flopy/mf6/data/mffileaccess.py | 4 +- flopy/utils/datafile.py | 2 +- 4 files changed, 193 insertions(+), 4 deletions(-) diff --git a/autotest/test_mf6.py b/autotest/test_mf6.py index 4aab5ccee..042df355c 100644 --- a/autotest/test_mf6.py +++ b/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 diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index a72d80524..d475bdd91 100644 --- a/autotest/test_model_splitter.py +++ b/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 diff --git a/flopy/mf6/data/mffileaccess.py b/flopy/mf6/data/mffileaccess.py index 2158b99bf..ce34a06f0 100644 --- a/flopy/mf6/data/mffileaccess.py +++ b/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)] diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 73a7f61b6..fb5906ae2 100644 --- a/flopy/utils/datafile.py +++ b/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":