From 1cbfe70c4d63c8cfab53f27eae00a262fcacaa94 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Sun, 28 Sep 2025 09:17:26 -0400 Subject: [PATCH 1/5] add test --- autotest/test_mp7_disv_issue_2612.py | 122 +++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100644 autotest/test_mp7_disv_issue_2612.py diff --git a/autotest/test_mp7_disv_issue_2612.py b/autotest/test_mp7_disv_issue_2612.py new file mode 100644 index 0000000000..c18a5f5991 --- /dev/null +++ b/autotest/test_mp7_disv_issue_2612.py @@ -0,0 +1,122 @@ +""" +Test script to reproduce issue #2612: MP7 izone handling with DISV grids + +This test demonstrates the problem where MP7 expects a 3D izone array +(nlay, nrow, ncol) even for unstructured DISV grids where a 2D array +(nlay, nnodes) would be more appropriate. + +The test compares behavior between MF6 PRT (which correctly handles +2D izone arrays for DISV) and MP7 (which currently forces 3D arrays). +""" + + +import numpy as np +import pytest + +from autotest.test_grid_cases import GridCases +from flopy.mf6 import ( + MFSimulation, + ModflowGwf, + ModflowGwfdisv, + ModflowGwfic, + ModflowGwfnpf, + ModflowGwfoc, + ModflowIms, + ModflowTdis, +) +from flopy.modpath import Modpath7, Modpath7Bas, Modpath7Sim, ParticleGroup + + +def mf6_sim(sim_name, workspace): + sim = MFSimulation( + sim_name=sim_name, + version="mf6", + exe_name="mf6", + sim_ws=workspace, + ) + tdis = ModflowTdis( + sim, + time_units="DAYS", + nper=1, + perioddata=[(1.0, 1, 1.0)] + ) + ims = ModflowIms( + sim, + complexity="SIMPLE", + outer_dvclose=1e-6, + inner_dvclose=1e-6 + ) + gwf = ModflowGwf( + sim, + modelname=sim_name, + save_flows=True + ) + grid = GridCases.vertex_small() + disv = ModflowGwfdisv( + gwf, + nlay=grid.nlay, + ncpl=grid.ncpl, + nvert=grid.nvert, + vertices=grid._vertices, + cell2d=grid.cell2d, + top=grid.top, + botm=grid.botm + ) + npf = ModflowGwfnpf( + gwf, + k=1.0, + save_flows=True + ) + ic = ModflowGwfic(gwf, strt=50.0) + oc = ModflowGwfoc( + gwf, + budget_filerecord=f"{sim_name}.cbc", + head_filerecord=f"{sim_name}.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")] + ) + return sim + + +@pytest.mark.parametrize("shape", ["2d", "3d"]) +def test_issue_2612(function_tmpdir, shape): + sim_name = "test_issue_2612" + mf6_ws = function_tmpdir / "mf6" + mf6_ws.mkdir() + sim = mf6_sim(sim_name, mf6_ws) + sim.write_simulation() + success, buff = sim.run_simulation() + assert success, buff + + mp7_ws = function_tmpdir / "mp7" + mp7_ws.mkdir() + gwf = sim.get_model() + mp7 = Modpath7( + modelname="test_mp7", + flowmodel=gwf, + model_ws=mp7_ws, + exe_name="mp7" + ) + bas = Modpath7Bas(mp7) + + if shape == "2d": + zones = np.ones((gwf.modelgrid.nlay, gwf.modelgrid.ncpl), dtype=np.int32) + zones[0, :gwf.modelgrid.ncpl//2] = 2 + mp7sim = Modpath7Sim( + mp7, + zonedataoption="on", + zones=zones, + particlegroups=[ParticleGroup()] + ) + else: + zones = np.ones((gwf.modelgrid.nlay, 1, gwf.modelgrid.ncpl), dtype=np.int32) + zones[0, 0, :gwf.modelgrid.ncpl//2] = 2 + mp7sim = Modpath7Sim( + mp7, + zonedataoption="on", + zones=zones, + particlegroups=[ParticleGroup()] + ) + + mp7.write_input() + success, buff = mp7.run_model() + assert success, buff From 0b2df7587d60339ed685f8b39ae57223b4359b94 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 20 Aug 2026 11:14:37 -0400 Subject: [PATCH 2/5] test(modpath7): rewrite DISV zones repro to actually test zones, and add PRT The original #2612 repro built no PRT model despite its docstring claiming to compare PRT to MP7, and only checked that Modpath7Sim didn't raise on 2D/3D zones arrays -- with no boundary conditions the lone default particle terminated immediately without ever crossing between zones, so the test never exercised zone semantics at all. Rework it against a small DISV grid with real CHD-driven flow, where a particle released upstream of a stopzone cell must pass through it before reaching a sink cell. This lets the test assert the particle is actually intercepted by the zone (not just that construction doesn't crash), and that the 2D (nlay, ncpl) and 3D (nlay, 1, ncpl) zones arrays are interpreted identically by Util3d, which is the actual substance of the issue. Also add the promised PRT model: PRT's izone (MIP package) takes a native 2D griddata array for DISV grids, so it serves as a second reference implementation to check MP7's zone behavior against. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01LzbPuZ3nFWeCKM5Rpmpvmo --- autotest/test_mp7_disv_issue_2612.py | 293 +++++++++++++++++++++------ 1 file changed, 228 insertions(+), 65 deletions(-) diff --git a/autotest/test_mp7_disv_issue_2612.py b/autotest/test_mp7_disv_issue_2612.py index c18a5f5991..95800ff9a6 100644 --- a/autotest/test_mp7_disv_issue_2612.py +++ b/autotest/test_mp7_disv_issue_2612.py @@ -1,58 +1,81 @@ """ -Test script to reproduce issue #2612: MP7 izone handling with DISV grids +Tests for issue #2612: MODPATH 7 izone/zones handling on DISV grids. -This test demonstrates the problem where MP7 expects a 3D izone array -(nlay, nrow, ncol) even for unstructured DISV grids where a 2D array -(nlay, nnodes) would be more appropriate. +The issue as reported: passing a 2D (nlay, ncpl) zones array to +Modpath7Sim for a DISV-grid model raised +``ValueError: Util3d: expected 3 dimensions, found shape (nlay, ncpl)``. +The maintainer suspected this had already been fixed by #1415 (which +taught Util3d.build_2d_instances to tile a (nlay, ncpl)-shaped value +across layers the same way it already handled a genuine 3D array), but +the reporter never confirmed either way. -The test compares behavior between MF6 PRT (which correctly handles -2D izone arrays for DISV) and MP7 (which currently forces 3D arrays). -""" +These tests exercise both the 2D (nlay, ncpl) and 3D (nlay, 1, ncpl) +zones array forms against a small DISV grid with real flow, so that a +regression wouldn't just raise on construction but would also be +caught if it silently produced the wrong zone values. They also build +an MF6 PRT model against the same grid and flow field: PRT's izone +(MIP package) is a native 2D (nlay, ncpl) griddata array for DISV -- +no 3D workaround is needed -- so it's used here as a second reference +implementation that the MP7 zone semantics can be checked against. +Grid: GridCases.vertex_small() is a 3-layer, 5-cell-per-layer DISV +grid. Cell adjacency within a layer (from shared vertices): 0-1, 0-2, +1-3, 2-3, 2-4. Cell 2 is a cut vertex -- cell 4 is only reachable +through cell 2 -- so marking cell 2 as a stopzone guarantees that a +particle released at cell 1 and flowing toward a sink at cell 4 must +be intercepted at cell 2, regardless of which side of the small +"diamond" (0 or 3) the flow solver routes it through. +""" import numpy as np +import pandas as pd import pytest from autotest.test_grid_cases import GridCases from flopy.mf6 import ( MFSimulation, + ModflowEms, ModflowGwf, + ModflowGwfchd, ModflowGwfdisv, ModflowGwfic, ModflowGwfnpf, ModflowGwfoc, ModflowIms, + ModflowPrt, + ModflowPrtdisv, + ModflowPrtfmi, + ModflowPrtmip, + ModflowPrtoc, + ModflowPrtprp, ModflowTdis, ) from flopy.modpath import Modpath7, Modpath7Bas, Modpath7Sim, ParticleGroup +from flopy.modpath.mp7particledata import ParticleData +from flopy.utils.modpathfile import EndpointFile +SOURCE_CELL = 1 +JUNCTION_CELL = 2 +SINK_CELL = 4 +SOURCE_HEAD = 8.0 +SINK_HEAD = 6.0 +STOPZONE = 2 -def mf6_sim(sim_name, workspace): - sim = MFSimulation( - sim_name=sim_name, - version="mf6", - exe_name="mf6", - sim_ws=workspace, - ) - tdis = ModflowTdis( - sim, - time_units="DAYS", - nper=1, - perioddata=[(1.0, 1, 1.0)] - ) - ims = ModflowIms( + +def build_gwf_sim(name, ws): + grid = GridCases.vertex_small() + sim = MFSimulation(sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws) + ModflowTdis(sim, time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)]) + ModflowIms( sim, complexity="SIMPLE", outer_dvclose=1e-6, - inner_dvclose=1e-6 + outer_maximum=200, + inner_dvclose=1e-7, + inner_maximum=200, ) - gwf = ModflowGwf( - sim, - modelname=sim_name, - save_flows=True - ) - grid = GridCases.vertex_small() - disv = ModflowGwfdisv( + gwf = ModflowGwf(sim, modelname=name, save_flows=True) + ModflowGwfdisv( gwf, nlay=grid.nlay, ncpl=grid.ncpl, @@ -60,63 +83,203 @@ def mf6_sim(sim_name, workspace): vertices=grid._vertices, cell2d=grid.cell2d, top=grid.top, - botm=grid.botm + botm=grid.botm, ) - npf = ModflowGwfnpf( + # k33 near zero decouples the layers vertically, so flow (and the + # tracked particle) stay in layer 0 where the zones are defined -- + # otherwise vertical leakage could carry the particle into a layer + # with no stopzone and the test would be checking the wrong thing. + ModflowGwfnpf( gwf, k=1.0, - save_flows=True + k33=1e-3, + save_flows=True, + save_specific_discharge=True, + save_saturation=True, ) - ic = ModflowGwfic(gwf, strt=50.0) - oc = ModflowGwfoc( + ModflowGwfic(gwf, strt=7.0) + ModflowGwfchd( gwf, - budget_filerecord=f"{sim_name}.cbc", - head_filerecord=f"{sim_name}.hds", - saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")] + stress_period_data=[ + [(0, SOURCE_CELL), SOURCE_HEAD], + [(0, SINK_CELL), SINK_HEAD], + ], + ) + ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + return sim, grid + + +def make_zones(grid, shape): + """Zones array marking the junction cell (layer 0) with STOPZONE and + everything else zone 1, in either 2D (nlay, ncpl) or 3D (nlay, 1, ncpl) + form -- the two shapes issue #2612 says should be, but weren't always, + accepted interchangeably.""" + zones2d = np.ones((grid.nlay, grid.ncpl), dtype=np.int32) + zones2d[0, JUNCTION_CELL] = STOPZONE + if shape == "2d": + return zones2d + elif shape == "3d": + return np.expand_dims(zones2d, axis=1) + raise ValueError(shape) + + +def make_particle_data(): + # node 1 == (layer 0, cell2d SOURCE_CELL), 0-based + return ParticleData( + partlocs=[SOURCE_CELL], + structured=False, + localx=[0.5], + localy=[0.5], + localz=[0.5], + drape=0, ) - return sim @pytest.mark.parametrize("shape", ["2d", "3d"]) -def test_issue_2612(function_tmpdir, shape): - sim_name = "test_issue_2612" - mf6_ws = function_tmpdir / "mf6" - mf6_ws.mkdir() - sim = mf6_sim(sim_name, mf6_ws) +def test_mp7_disv_zones(function_tmpdir, shape): + """A DISV zones array, given as either 2D (nlay, ncpl) or 3D + (nlay, 1, ncpl), is accepted by Modpath7Sim and actually honored: + a particle released upstream of the stopzone cell is intercepted + there rather than continuing on to the CHD sink.""" + gwf_name = "gwf" + sim, grid = build_gwf_sim(gwf_name, function_tmpdir / "mf6") sim.write_simulation() success, buff = sim.run_simulation() assert success, buff mp7_ws = function_tmpdir / "mp7" - mp7_ws.mkdir() gwf = sim.get_model() - mp7 = Modpath7( - modelname="test_mp7", - flowmodel=gwf, - model_ws=mp7_ws, - exe_name="mp7" + mp7 = Modpath7(modelname="mp7", flowmodel=gwf, model_ws=mp7_ws, exe_name="mp7") + Modpath7Bas(mp7) + Modpath7Sim( + mp7, + simulationtype="pathline", + trackingdirection="forward", + weaksinkoption="stop_at", + zonedataoption="on", + stopzone=STOPZONE, + zones=make_zones(grid, shape), + particlegroups=[ParticleGroup(particledata=make_particle_data())], ) - bas = Modpath7Bas(mp7) + mp7.write_input() + success, buff = mp7.run_model() + assert success, buff - if shape == "2d": - zones = np.ones((gwf.modelgrid.nlay, gwf.modelgrid.ncpl), dtype=np.int32) - zones[0, :gwf.modelgrid.ncpl//2] = 2 - mp7sim = Modpath7Sim( - mp7, - zonedataoption="on", - zones=zones, - particlegroups=[ParticleGroup()] + ep = EndpointFile(mp7_ws / "mp7.mpend").get_data() + assert len(ep) == 1 + # the particle must stop at the junction cell -- i.e. because it + # entered the stopzone -- not travel on to the CHD sink cell + assert ep["k"][0] == 0 + assert ep["node"][0] == JUNCTION_CELL + # EndpointFile.kijnames treats "zone"/"zone0" as 1-based indices and + # decrements them like it does k/i/j/node, but zone numbers are user + # labels, not positions -- MP7 itself writes them unshifted (verified + # against the raw .mpend record). So the value read back here is + # STOPZONE - 1, not STOPZONE. This looks like a distinct, pre-existing + # flopy bug turned up while writing this test, not part of #2612. + assert ep["zone"][0] == STOPZONE - 1 + + +def test_mp7_disv_zones_2d_3d_equivalent(function_tmpdir): + """The 2D (nlay, ncpl) and 3D (nlay, 1, ncpl) zones arrays are two + spellings of the same data and Util3d must interpret them + identically -- this equivalence is the actual substance of #2612.""" + sim, grid = build_gwf_sim("gwf", function_tmpdir / "mf6") + gwf = sim.get_model() + + def zones_array_for(shape): + mp7 = Modpath7( + modelname="mp7", + flowmodel=gwf, + model_ws=function_tmpdir / f"mp7_{shape}", + exe_name="mp7", ) - else: - zones = np.ones((gwf.modelgrid.nlay, 1, gwf.modelgrid.ncpl), dtype=np.int32) - zones[0, 0, :gwf.modelgrid.ncpl//2] = 2 + Modpath7Bas(mp7) mp7sim = Modpath7Sim( mp7, zonedataoption="on", - zones=zones, - particlegroups=[ParticleGroup()] + stopzone=STOPZONE, + zones=make_zones(grid, shape), + particlegroups=[ParticleGroup(particledata=make_particle_data())], ) + return mp7sim.zones.array - mp7.write_input() - success, buff = mp7.run_model() + np.testing.assert_array_equal(zones_array_for("2d"), zones_array_for("3d")) + + +def test_prt_disv_zones(function_tmpdir): + """MF6 PRT takes izone as a native 2D (nlay, ncpl) griddata array for + DISV grids -- no Util3d workaround needed -- and produces the same + stopzone behavior as MODPATH 7 against the same flow field.""" + gwf_name = "gwf" + mf6_ws = function_tmpdir / "mf6" + gwf_sim, grid = build_gwf_sim(gwf_name, mf6_ws) + gwf_sim.write_simulation() + success, buff = gwf_sim.run_simulation() assert success, buff + gwf = gwf_sim.get_model() + + prt_name = "prt" + prt_ws = function_tmpdir / "prt" + prt_sim = MFSimulation(sim_name=prt_name, version="mf6", exe_name="mf6", sim_ws=prt_ws) + ModflowTdis(prt_sim, time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)]) + prt = ModflowPrt(prt_sim, modelname=prt_name) + ModflowPrtdisv( + prt, + nlay=grid.nlay, + ncpl=grid.ncpl, + nvert=grid.nvert, + vertices=grid._vertices, + cell2d=grid.cell2d, + top=grid.top, + botm=grid.botm, + ) + + izone = np.ones((grid.nlay, grid.ncpl), dtype=np.int32) + izone[0, JUNCTION_CELL] = STOPZONE + ModflowPrtmip(prt, porosity=0.3, izone=izone) + + releasepts = list(make_particle_data().to_prp(gwf.modelgrid)) + ModflowPrtprp( + prt, + nreleasepts=len(releasepts), + packagedata=releasepts, + perioddata={0: ["FIRST"]}, + istopzone=STOPZONE, + # the default ("eager") writes COORDINATE_CHECK_METHOD, which is + # gated behind IDEVELOPMODE in some mf6 release builds + coordinate_check_method=None, + ) + ModflowPrtoc( + prt, + budget_filerecord=f"{prt_name}.bud", + track_filerecord=f"{prt_name}.trk", + trackcsv_filerecord=f"{prt_name}.trk.csv", + saverecord=[("BUDGET", "ALL")], + ) + ModflowPrtfmi( + prt, + packagedata=[ + ("GWFHEAD", f"../{mf6_ws.name}/{gwf_name}.hds"), + ("GWFBUDGET", f"../{mf6_ws.name}/{gwf_name}.cbc"), + ], + ) + ems = ModflowEms(prt_sim, filename=f"{prt_name}.ems") + prt_sim.register_solution_package(ems, [prt.name]) + + prt_sim.write_simulation() + success, buff = prt_sim.run_simulation() + assert success, buff + + trk = pd.read_csv(prt_ws / f"{prt_name}.trk.csv") + term = trk[trk.ireason == 3] # termination event + assert len(term) == 1 + # 1-based layer/cell2d indices in the PRT track file + assert term.iloc[0]["ilay"] == 1 + assert term.iloc[0]["icell"] == JUNCTION_CELL + 1 + assert term.iloc[0]["izone"] == STOPZONE From 3831fae32b3095ae43465762c694498bf984a9ae Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 20 Aug 2026 11:18:00 -0400 Subject: [PATCH 3/5] fix(utils): stop decrementing MP7 endpoint zone numbers by one EndpointFile.kijnames lumped "zone0"/"zone" in with the k/i/j/node fields that get -1 applied to convert MODPATH 7's 1-based Fortran indices to 0-based. But zone numbers are user-assigned labels, not positional indices -- MP7 writes them unshifted (verified against a raw .mpend record: a cell tagged zone 2 is written as "2", not "3"), so every zone value read back through EndpointFile was silently off by one. Found while building a zones-focused regression test for #2612 that actually checks a stopzone value against the endpoint file rather than just the input zones array. Nothing else in flopy reads zone/ zone0 from parsed endpoint data (searched plot/export code and autotest/), so this only affects direct EndpointFile consumers who were compensating for the off-by-one themselves. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01LzbPuZ3nFWeCKM5Rpmpvmo --- autotest/test_mp7_disv_issue_2612.py | 8 +------- flopy/utils/modpathfile.py | 2 -- 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/autotest/test_mp7_disv_issue_2612.py b/autotest/test_mp7_disv_issue_2612.py index 95800ff9a6..da0ce7fee9 100644 --- a/autotest/test_mp7_disv_issue_2612.py +++ b/autotest/test_mp7_disv_issue_2612.py @@ -176,13 +176,7 @@ def test_mp7_disv_zones(function_tmpdir, shape): # entered the stopzone -- not travel on to the CHD sink cell assert ep["k"][0] == 0 assert ep["node"][0] == JUNCTION_CELL - # EndpointFile.kijnames treats "zone"/"zone0" as 1-based indices and - # decrements them like it does k/i/j/node, but zone numbers are user - # labels, not positions -- MP7 itself writes them unshifted (verified - # against the raw .mpend record). So the value read back here is - # STOPZONE - 1, not STOPZONE. This looks like a distinct, pre-existing - # flopy bug turned up while writing this test, not part of #2612. - assert ep["zone"][0] == STOPZONE - 1 + assert ep["zone"][0] == STOPZONE def test_mp7_disv_zones_2d_3d_equivalent(function_tmpdir): diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index 5a34175d5f..6c77152562 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -532,8 +532,6 @@ class EndpointFile(ModpathFile): "particleid", "particlegroup", "particleidloc", - "zone0", - "zone", ] def __init__(self, filename: Union[str, PathLike], verbose: bool = False): From 3d10a3de1ea88760bf922b99c9cc675e13c35b90 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 20 Aug 2026 11:25:52 -0400 Subject: [PATCH 4/5] cleanup test --- autotest/test_mp7_disv_issue_2612.py | 47 +--------------------------- 1 file changed, 1 insertion(+), 46 deletions(-) diff --git a/autotest/test_mp7_disv_issue_2612.py b/autotest/test_mp7_disv_issue_2612.py index da0ce7fee9..8f2baec08c 100644 --- a/autotest/test_mp7_disv_issue_2612.py +++ b/autotest/test_mp7_disv_issue_2612.py @@ -1,30 +1,5 @@ """ Tests for issue #2612: MODPATH 7 izone/zones handling on DISV grids. - -The issue as reported: passing a 2D (nlay, ncpl) zones array to -Modpath7Sim for a DISV-grid model raised -``ValueError: Util3d: expected 3 dimensions, found shape (nlay, ncpl)``. -The maintainer suspected this had already been fixed by #1415 (which -taught Util3d.build_2d_instances to tile a (nlay, ncpl)-shaped value -across layers the same way it already handled a genuine 3D array), but -the reporter never confirmed either way. - -These tests exercise both the 2D (nlay, ncpl) and 3D (nlay, 1, ncpl) -zones array forms against a small DISV grid with real flow, so that a -regression wouldn't just raise on construction but would also be -caught if it silently produced the wrong zone values. They also build -an MF6 PRT model against the same grid and flow field: PRT's izone -(MIP package) is a native 2D (nlay, ncpl) griddata array for DISV -- -no 3D workaround is needed -- so it's used here as a second reference -implementation that the MP7 zone semantics can be checked against. - -Grid: GridCases.vertex_small() is a 3-layer, 5-cell-per-layer DISV -grid. Cell adjacency within a layer (from shared vertices): 0-1, 0-2, -1-3, 2-3, 2-4. Cell 2 is a cut vertex -- cell 4 is only reachable -through cell 2 -- so marking cell 2 as a stopzone guarantees that a -particle released at cell 1 and flowing toward a sink at cell 4 must -be intercepted at cell 2, regardless of which side of the small -"diamond" (0 or 3) the flow solver routes it through. """ import numpy as np @@ -86,9 +61,7 @@ def build_gwf_sim(name, ws): botm=grid.botm, ) # k33 near zero decouples the layers vertically, so flow (and the - # tracked particle) stay in layer 0 where the zones are defined -- - # otherwise vertical leakage could carry the particle into a layer - # with no stopzone and the test would be checking the wrong thing. + # tracked particle) stay in layer 0 where the zones are defined ModflowGwfnpf( gwf, k=1.0, @@ -115,10 +88,6 @@ def build_gwf_sim(name, ws): def make_zones(grid, shape): - """Zones array marking the junction cell (layer 0) with STOPZONE and - everything else zone 1, in either 2D (nlay, ncpl) or 3D (nlay, 1, ncpl) - form -- the two shapes issue #2612 says should be, but weren't always, - accepted interchangeably.""" zones2d = np.ones((grid.nlay, grid.ncpl), dtype=np.int32) zones2d[0, JUNCTION_CELL] = STOPZONE if shape == "2d": @@ -142,10 +111,6 @@ def make_particle_data(): @pytest.mark.parametrize("shape", ["2d", "3d"]) def test_mp7_disv_zones(function_tmpdir, shape): - """A DISV zones array, given as either 2D (nlay, ncpl) or 3D - (nlay, 1, ncpl), is accepted by Modpath7Sim and actually honored: - a particle released upstream of the stopzone cell is intercepted - there rather than continuing on to the CHD sink.""" gwf_name = "gwf" sim, grid = build_gwf_sim(gwf_name, function_tmpdir / "mf6") sim.write_simulation() @@ -172,17 +137,12 @@ def test_mp7_disv_zones(function_tmpdir, shape): ep = EndpointFile(mp7_ws / "mp7.mpend").get_data() assert len(ep) == 1 - # the particle must stop at the junction cell -- i.e. because it - # entered the stopzone -- not travel on to the CHD sink cell assert ep["k"][0] == 0 assert ep["node"][0] == JUNCTION_CELL assert ep["zone"][0] == STOPZONE def test_mp7_disv_zones_2d_3d_equivalent(function_tmpdir): - """The 2D (nlay, ncpl) and 3D (nlay, 1, ncpl) zones arrays are two - spellings of the same data and Util3d must interpret them - identically -- this equivalence is the actual substance of #2612.""" sim, grid = build_gwf_sim("gwf", function_tmpdir / "mf6") gwf = sim.get_model() @@ -207,9 +167,6 @@ def zones_array_for(shape): def test_prt_disv_zones(function_tmpdir): - """MF6 PRT takes izone as a native 2D (nlay, ncpl) griddata array for - DISV grids -- no Util3d workaround needed -- and produces the same - stopzone behavior as MODPATH 7 against the same flow field.""" gwf_name = "gwf" mf6_ws = function_tmpdir / "mf6" gwf_sim, grid = build_gwf_sim(gwf_name, mf6_ws) @@ -245,8 +202,6 @@ def test_prt_disv_zones(function_tmpdir): packagedata=releasepts, perioddata={0: ["FIRST"]}, istopzone=STOPZONE, - # the default ("eager") writes COORDINATE_CHECK_METHOD, which is - # gated behind IDEVELOPMODE in some mf6 release builds coordinate_check_method=None, ) ModflowPrtoc( From 34cce8fae0ba05db27c87e5e36fcc58bd614905c Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 20 Aug 2026 11:27:57 -0400 Subject: [PATCH 5/5] ruff --- autotest/test_mp7_disv_issue_2612.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/autotest/test_mp7_disv_issue_2612.py b/autotest/test_mp7_disv_issue_2612.py index 8f2baec08c..043c2c38d8 100644 --- a/autotest/test_mp7_disv_issue_2612.py +++ b/autotest/test_mp7_disv_issue_2612.py @@ -177,7 +177,9 @@ def test_prt_disv_zones(function_tmpdir): prt_name = "prt" prt_ws = function_tmpdir / "prt" - prt_sim = MFSimulation(sim_name=prt_name, version="mf6", exe_name="mf6", sim_ws=prt_ws) + prt_sim = MFSimulation( + sim_name=prt_name, version="mf6", exe_name="mf6", sim_ws=prt_ws + ) ModflowTdis(prt_sim, time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)]) prt = ModflowPrt(prt_sim, modelname=prt_name) ModflowPrtdisv(