From da25429a846275a1408d9afe0967d9d099a7147e Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Sun, 23 Aug 2026 09:38:05 -0500 Subject: [PATCH 1/2] test(model_splitter): cover the horizontal flow barrier package The HFB package had no splitting test. Split a model with barriers on either side of the split column and check that each model keeps its own barriers and that the heads reconstruct to the original model, then check that splitting a model along a barrier is an error. --- autotest/test_model_splitter.py | 69 +++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 6bd90dff8..ab09f61c4 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -174,6 +174,75 @@ 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_hfb_model_splitter(function_tmpdir): + nlay, nrow, ncol = 1, 10, 10 + + sim = flopy.mf6.MFSimulation(sim_name="hfb", sim_ws=function_tmpdir, exe_name="mf6") + flopy.mf6.ModflowTdis(sim) + flopy.mf6.ModflowIms(sim, complexity="simple") + gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb", save_flows=True) + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=100.0, + delc=100.0, + top=10.0, + botm=[0.0], + ) + flopy.mf6.ModflowGwfnpf(gwf, k=1.0) + flopy.mf6.ModflowGwfic(gwf, strt=10.0) + chd = [[(0, i, 0), 10.0] for i in range(nrow)] + chd += [[(0, i, ncol - 1), 1.0] for i in range(nrow)] + flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd) + + # a barrier on either side of the column the model is split on + hfb = [[(0, i, 2), (0, i, 3), 1e-5] for i in range(nrow)] + hfb += [[(0, i, 6), (0, i, 7), 1e-5] for i in range(nrow)] + flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfb}) + flopy.mf6.ModflowGwfoc(gwf, head_filerecord="hfb.hds", saverecord=[("HEAD", "ALL")]) + sim.write_simulation() + sim.run_simulation() + + original_heads = gwf.output.head().get_alldata()[-1] + + array = np.zeros((nrow, ncol), dtype=int) + array[:, 5:] = 1 + + mfsplit = Mf6Splitter(sim) + new_sim = mfsplit.split_model(array) + new_sim.set_sim_path(function_tmpdir / "split_model") + new_sim.write_simulation() + new_sim.run_simulation() + + heads = {} + for mkey in (0, 1): + ml = new_sim.get_model(f"hfb_{mkey}") + spd = ml.get_package("hfb").stress_period_data.get_data(0) + if len(spd) != nrow: + raise AssertionError( + f"Model {mkey} has {len(spd)} barriers, expected {nrow}" + ) + + heads[mkey] = ml.output.head().get_alldata()[-1] + + new_heads = mfsplit.reconstruct_array(heads) + + err_msg = "Heads from original and split models do not match" + np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg) + + # a barrier cannot be split, both of its cells must be in one model + sim = MFSimulation.load(sim_ws=function_tmpdir) + array = np.zeros((nrow, ncol), dtype=int) + array[:, 3:] = 1 + + mfsplit = Mf6Splitter(sim) + with pytest.raises(AssertionError, match="split along faults"): + mfsplit.split_model(array) + + @requires_exe("mf6") @requires_pkg("pymetis") @pytest.mark.slow From e7cc42d5a85216673cbd2b3f6290ea8ff38ef565 Mon Sep 17 00:00:00 2001 From: Joseph Hughes Date: Sun, 23 Aug 2026 09:55:20 -0500 Subject: [PATCH 2/2] test(model_splitter): add a barrier junction and vertical barriers Cover the corner four cells share, where two barriers between east west neighbors and two between north south neighbors meet, and check that splitting on either the row or the column the junction straddles is an error. Add barriers on vertical connections, which the splitter never divides because both cells of one are in the same stack. --- autotest/test_model_splitter.py | 73 ++++++++++++++++++++++----------- 1 file changed, 48 insertions(+), 25 deletions(-) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index ab09f61c4..1872d9f21 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -176,11 +176,13 @@ def test_model_with_lak_sfr_mvr(function_tmpdir): @requires_exe("mf6") def test_hfb_model_splitter(function_tmpdir): - nlay, nrow, ncol = 1, 10, 10 + nlay, nrow, ncol = 3, 10, 10 sim = flopy.mf6.MFSimulation(sim_name="hfb", sim_ws=function_tmpdir, exe_name="mf6") flopy.mf6.ModflowTdis(sim) - flopy.mf6.ModflowIms(sim, complexity="simple") + flopy.mf6.ModflowIms( + sim, complexity="simple", outer_dvclose=1e-9, inner_dvclose=1e-10 + ) gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb", save_flows=True) flopy.mf6.ModflowGwfdis( gwf, @@ -189,18 +191,27 @@ def test_hfb_model_splitter(function_tmpdir): ncol=ncol, delr=100.0, delc=100.0, - top=10.0, - botm=[0.0], + top=30.0, + botm=[20.0, 10.0, 0.0], ) flopy.mf6.ModflowGwfnpf(gwf, k=1.0) - flopy.mf6.ModflowGwfic(gwf, strt=10.0) - chd = [[(0, i, 0), 10.0] for i in range(nrow)] - chd += [[(0, i, ncol - 1), 1.0] for i in range(nrow)] + flopy.mf6.ModflowGwfic(gwf, strt=25.0) + chd = [[(0, i, 0), 25.0] for i in range(nrow)] + chd += [[(0, i, ncol - 1), 20.0] for i in range(nrow)] flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd) - # a barrier on either side of the column the model is split on - hfb = [[(0, i, 2), (0, i, 3), 1e-5] for i in range(nrow)] - hfb += [[(0, i, 6), (0, i, 7), 1e-5] for i in range(nrow)] + # four barriers forming a plus on the corner shared by the cells + # (4, 4), (4, 5), (5, 4), and (5, 5): two between east west neighbors + # and two between north south neighbors + hfb = [ + [(0, 4, 4), (0, 4, 5), 1e-5], # east west, north of the corner + [(0, 5, 4), (0, 5, 5), 1e-5], # east west, south of the corner + [(0, 4, 4), (0, 5, 4), 1e-5], # north south, west of the corner + [(0, 4, 5), (0, 5, 5), 1e-5], # north south, east of the corner + ] + # a barrier on a vertical connection, which the splitter never divides + # because both of its cells are in the same stack + hfb += [[(0, i, j), (1, i, j), 1e-5] for i in range(3) for j in range(7, 10)] flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfb}) flopy.mf6.ModflowGwfoc(gwf, head_filerecord="hfb.hds", saverecord=[("HEAD", "ALL")]) sim.write_simulation() @@ -208,8 +219,9 @@ def test_hfb_model_splitter(function_tmpdir): original_heads = gwf.output.head().get_alldata()[-1] + # split away from every barrier array = np.zeros((nrow, ncol), dtype=int) - array[:, 5:] = 1 + array[:, 7:] = 1 mfsplit = Mf6Splitter(sim) new_sim = mfsplit.split_model(array) @@ -218,29 +230,40 @@ def test_hfb_model_splitter(function_tmpdir): new_sim.run_simulation() heads = {} + nbarrier = 0 for mkey in (0, 1): ml = new_sim.get_model(f"hfb_{mkey}") - spd = ml.get_package("hfb").stress_period_data.get_data(0) - if len(spd) != nrow: - raise AssertionError( - f"Model {mkey} has {len(spd)} barriers, expected {nrow}" - ) + pkg = ml.get_package("hfb") + if pkg is None: + raise AssertionError(f"Model {mkey} has no HFB package") + nbarrier += len(pkg.stress_period_data.get_data(0)) heads[mkey] = ml.output.head().get_alldata()[-1] + if nbarrier != len(hfb): + raise AssertionError( + f"Split models have {nbarrier} barriers, expected {len(hfb)}" + ) + new_heads = mfsplit.reconstruct_array(heads) err_msg = "Heads from original and split models do not match" - np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg) - - # a barrier cannot be split, both of its cells must be in one model - sim = MFSimulation.load(sim_ws=function_tmpdir) - array = np.zeros((nrow, ncol), dtype=int) - array[:, 3:] = 1 + np.testing.assert_allclose(new_heads, original_heads, atol=1e-6, err_msg=err_msg) + + # a barrier cannot be split, both of its cells must be in one model. + # splitting on the column the plus straddles cuts its east west barriers + # and splitting on the row cuts its north south barriers + for axis in (0, 1): + sim = MFSimulation.load(sim_ws=function_tmpdir) + array = np.zeros((nrow, ncol), dtype=int) + if axis == 0: + array[5:, :] = 1 + else: + array[:, 5:] = 1 - mfsplit = Mf6Splitter(sim) - with pytest.raises(AssertionError, match="split along faults"): - mfsplit.split_model(array) + mfsplit = Mf6Splitter(sim) + with pytest.raises(AssertionError, match="split along faults"): + mfsplit.split_model(array) @requires_exe("mf6")