From 5bc119813c3147a28dbbbfac3bc44b205ca67f66 Mon Sep 17 00:00:00 2001 From: Bonelli Date: Tue, 24 Feb 2026 09:46:04 -0500 Subject: [PATCH 01/10] feat: support nwt to mf6 binary output file conversion --- autotest/test_nwt_to_mf6.py | 668 +++++++++++++++++++++++++++++ autotest/test_postprocessing.py | 255 +++++++++++ flopy/mf6/utils/__init__.py | 7 +- flopy/mf6/utils/binarygrid_util.py | 310 +++++++++++++ flopy/mf6/utils/postprocessing.py | 272 +++++++++++- flopy/utils/__init__.py | 1 + flopy/utils/nwt_to_mf6.py | 521 ++++++++++++++++++++++ flopy/utils/postprocessing.py | 197 +++++++++ 8 files changed, 2211 insertions(+), 20 deletions(-) create mode 100644 autotest/test_nwt_to_mf6.py create mode 100644 flopy/utils/nwt_to_mf6.py diff --git a/autotest/test_nwt_to_mf6.py b/autotest/test_nwt_to_mf6.py new file mode 100644 index 0000000000..7376f1f2ee --- /dev/null +++ b/autotest/test_nwt_to_mf6.py @@ -0,0 +1,668 @@ +import numpy as np +import pytest +from modflow_devtools.markers import requires_exe + + +def test_get_icelltype_from_laytyp(): + from flopy.utils.nwt_to_mf6 import get_icelltype_from_laytyp + + # 1D array + laytyp = np.array([1, 0, 0, 1]) + icelltype = get_icelltype_from_laytyp(laytyp) + assert icelltype.shape == laytyp.shape + assert np.array_equal(icelltype, [1, 0, 0, 1]) + + # scalar + laytyp = 0 + icelltype = get_icelltype_from_laytyp(laytyp) + assert icelltype == 0 + + laytyp = 1 + icelltype = get_icelltype_from_laytyp(laytyp) + assert icelltype == 1 + + # various LAYTYP values + laytyp = np.array([0, 1, 2, 3, -1]) + icelltype = get_icelltype_from_laytyp(laytyp) + # All > 0 map to 1, 0 stays 0, <0 stays 0 + assert np.array_equal(icelltype, [0, 1, 1, 1, 0]) + + +def test_mfgrdfile_write_roundtrip(tmp_path): + from flopy.mf6.utils import MfGrdFile, get_structured_connectivity + + nlay, nrow, ncol = 2, 3, 4 + delr = np.ones(ncol) * 100.0 + delc = np.ones(nrow) * 50.0 + top = np.ones((nrow, ncol)) * 10.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 5.0 + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) + icelltype[0] = 1 # top layer convertible + + grb_file = tmp_path / "test.dis.grb" + MfGrdFile.write_dis( + str(grb_file), + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + ia, + ja, + icelltype=icelltype, + ) + + grb = MfGrdFile(str(grb_file)) + + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + assert grb.nodes == nlay * nrow * ncol + assert grb.nja == nja + + np.testing.assert_array_almost_equal(grb.delr, delr) + np.testing.assert_array_almost_equal(grb.delc, delc) + + # Build expected TOP for all cells (MF6 format) + # Layer 1: model top, Layer 2+: bottom of layer above + # In Fortran order (layer-interleaved): [L0[0,0], L1[0,0], L2[0,0], L0[0,1], ...] + top_expected = np.zeros(nlay * nrow * ncol) + top_flat = top.flatten(order="F") + botm_flat = botm.flatten(order="F") + + for i in range(nrow * ncol): + # Layer 0: use model top + top_expected[i * nlay] = top_flat[i] + # Layers 1+: use bottom of layer above + for k in range(1, nlay): + top_expected[i * nlay + k] = botm_flat[i * nlay + (k - 1)] + + np.testing.assert_array_almost_equal(grb.top, top_expected) + np.testing.assert_array_almost_equal(grb.bot, botm.flatten(order="F")) + + np.testing.assert_array_equal(grb.ia, ia) + np.testing.assert_array_equal(grb.ja, ja) + + icelltype_read = grb._datadict["ICELLTYPE"] + np.testing.assert_array_equal(icelltype_read, icelltype.flatten(order="F")) + + +def test_mfgrdfile_write_with_idomain(tmp_path): + from flopy.mf6.utils import MfGrdFile, get_structured_connectivity + + nlay, nrow, ncol = 1, 3, 3 + delr = np.ones(ncol) + delc = np.ones(nrow) + top = np.ones((nrow, ncol)) + botm = np.zeros((nlay, nrow, ncol)) + + # center cell inactive + idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) + idomain[0, 1, 1] = 0 + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol, idomain) + + grb_file = tmp_path / "test_idomain.dis.grb" + MfGrdFile.write_dis( + str(grb_file), nlay, nrow, ncol, delr, delc, top, botm, ia, ja, idomain=idomain + ) + + grb = MfGrdFile(str(grb_file)) + idomain_read = grb._datadict["IDOMAIN"] + np.testing.assert_array_equal(idomain_read, idomain.flatten(order="F")) + + +def test_mfgrdfile_write_validation(): + from flopy.mf6.utils import MfGrdFile, get_structured_connectivity + + nlay, nrow, ncol = 1, 2, 2 + delr = np.ones(ncol) + delc = np.ones(nrow) + top = np.ones((nrow, ncol)) + botm = np.zeros((nlay, nrow, ncol)) + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + with pytest.raises(ValueError, match="delr length"): + MfGrdFile.write_dis( + "test.grb", + nlay, + nrow, + ncol, + np.ones(3), # Wrong length + delc, + top, + botm, + ia, + ja, + ) + + with pytest.raises(ValueError, match="ia length"): + MfGrdFile.write_dis( + "test.grb", + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + np.ones(10), # Wrong length (should be ncells + 1 = 5) + ja, + ) + + +def test_nwt_to_mf6_converter_init(function_tmpdir): + import numpy as np + + from flopy.utils import NwtToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + nlay, nrow, ncol = 2, 3, 4 + delr = np.ones(ncol) * 100.0 + delc = np.ones(nrow) * 100.0 + top = np.ones((nrow, ncol)) * 10.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 5.0 + laytyp = np.array([1, 0]) # Top layer convertible + + hds_file = function_tmpdir / "test.hds" + cbc_file = function_tmpdir / "test.cbc" + + head_data = np.ones((nlay, nrow, ncol)) * 8.0 + HeadFile.write( + str(hds_file), + {(1, 1): head_data}, + nlay=nlay, + nrow=nrow, + ncol=ncol, + ) + + from flopy.mf6.utils import get_structured_connectivity, get_structured_faceflows + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + flowja = np.ones(nja) * 0.5 + + frf, fff, flf = get_structured_faceflows( + flowja, grb_file=None, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) + + bud_data = [ + { + "data": frf.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW RIGHT FACE", + "imeth": 1, + }, + { + "data": fff.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW FRONT FACE", + "imeth": 1, + }, + { + "data": flf.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW LOWER FACE", + "imeth": 1, + }, + ] + + CellBudgetFile.write( + str(cbc_file), + bud_data, + nlay=nlay, + nrow=nrow, + ncol=ncol, + ) + + converter = NwtToMf6Converter( + str(hds_file), + str(cbc_file), + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + laytyp, + model_ws=function_tmpdir, + ) + + assert converter.nlay == nlay + assert converter.nrow == nrow + assert converter.ncol == ncol + assert converter.ncells == nlay * nrow * ncol + assert len(converter.times) == 1 + assert len(converter.kstpkper) == 1 + + assert converter.icelltype.shape == (nlay,) + assert np.array_equal(converter.icelltype, [1, 0]) + assert converter.icelltype_3d.shape == (nlay, nrow, ncol) + assert np.all(converter.icelltype_3d[0] == 1) + assert np.all(converter.icelltype_3d[1] == 0) + + +def test_nwt_to_mf6_converter_convert(function_tmpdir): + import numpy as np + + from flopy.mf6.utils import ( + MfGrdFile, + get_structured_connectivity, + get_structured_faceflows, + ) + from flopy.utils import NwtToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + nlay, nrow, ncol = 2, 3, 4 + delr = np.ones(ncol) * 100.0 + delc = np.ones(nrow) * 100.0 + top = np.ones((nrow, ncol)) * 10.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 5.0 + laytyp = np.array([1, 0]) + + hds_file = function_tmpdir / "test.hds" + cbc_file = function_tmpdir / "test.cbc" + + head_data = np.ones((nlay, nrow, ncol)) * 8.0 + HeadFile.write(str(hds_file), {(1, 1): head_data}, nlay=nlay, nrow=nrow, ncol=ncol) + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + flowja = np.ones(nja) * 0.5 + frf, fff, flf = get_structured_faceflows( + flowja, grb_file=None, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) + + bud_data = [ + { + "data": frf.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW RIGHT FACE", + "imeth": 1, + }, + { + "data": fff.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW FRONT FACE", + "imeth": 1, + }, + { + "data": flf.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW LOWER FACE", + "imeth": 1, + }, + ] + CellBudgetFile.write(str(cbc_file), bud_data, nlay=nlay, nrow=nrow, ncol=ncol) + + converter = NwtToMf6Converter( + str(hds_file), + str(cbc_file), + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + laytyp, + model_ws=function_tmpdir, + ) + + output_dir = function_tmpdir / "mf6_output" + result = converter.convert(str(output_dir), verbose=False) + + assert result["grb"].exists() + assert result["hds"].exists() + assert result["bud"].exists() + + grb = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + hds_mf6 = HeadFile(str(result["hds"])) + head_read = hds_mf6.get_data(idx=0) # Read first record + np.testing.assert_array_almost_equal(head_read, head_data) + + bud_mf6 = CellBudgetFile(str(result["bud"])) + texts = bud_mf6.textlist + # Convert to strings and strip for comparison + texts_str = [ + t.decode().strip() if isinstance(t, bytes) else str(t).strip() for t in texts + ] + assert "FLOW-JA-FACE" in texts_str + assert "DATA-SAT" in texts_str + # DATA-SPDIS skipped for now + + +@requires_exe("mfnwt") +@pytest.mark.slow +def test_nwt_to_mf6_watertable_model(function_tmpdir): + from flopy.mf6.utils import MfGrdFile + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowGhb, + ModflowNwt, + ModflowOc, + ModflowRch, + ModflowUpw, + ) + from flopy.utils import NwtToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + modelname = "watertable" + + nlay, nrow, ncol = 1, 1, 100 + delr = 50.0 + delc = 1.0 + + h1, h2 = 20.0, 11.0 + + top = 25.0 + botm = 0.0 + hk = 50.0 + + strt = np.zeros((nlay, nrow, ncol), dtype=float) + strt[0, 0, 0] = h1 + strt[0, 0, -1] = h2 + + rchrate = 0.001 + + h_adj1 = h1 - (h1 - h2) / ncol + h_adj2 = h2 + (h1 - h2) / ncol + + b1 = 0.5 * (h1 + h_adj1) + b2 = 0.5 * (h2 + h_adj2) + c1 = hk * b1 * delc / (0.5 * delr) + c2 = hk * b2 * delc / (0.5 * delr) + + ghb_dtype = ModflowGhb.get_default_dtype() + stress_period_data = np.zeros((2), dtype=ghb_dtype) + stress_period_data = stress_period_data.view(np.recarray) + stress_period_data[0] = (0, 0, 0, h1, c1) + stress_period_data[1] = (0, 0, ncol - 1, h2, c2) + + mf = Modflow( + modelname=modelname, + exe_name="mfnwt", + model_ws=function_tmpdir, + version="mfnwt", + ) + ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + perlen=1, + nstp=1, + steady=True, + ) + ModflowBas(mf, ibound=1, strt=strt) + ModflowUpw(mf, hk=hk, laytyp=1, ipakcb=53) # laytyp=1 for convertible + ModflowGhb(mf, stress_period_data=stress_period_data) + ModflowRch(mf, rech=rchrate, nrchop=1) + oc = ModflowOc( + mf, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, + ) + oc.reset_budgetunit(budgetunit=53, fname=f"{modelname}.cbc") + ModflowNwt(mf) + + mf.write_input() + success, _ = mf.run_model(silent=True) + assert success, "NWT model run failed" + + hds_file = function_tmpdir / f"{modelname}.hds" + cbc_file = function_tmpdir / f"{modelname}.cbc" + assert hds_file.exists(), "Head file not created" + assert cbc_file.exists(), "Budget file not created" + + converter = NwtToMf6Converter( + str(hds_file), + str(cbc_file), + nlay, + nrow, + ncol, + delr=np.ones(ncol) * delr, + delc=np.ones(nrow) * delc, + top=np.ones((nrow, ncol)) * top, + botm=np.ones((nlay, nrow, ncol)) * botm, + laytyp=np.array([1]), # Convertible + model_ws=function_tmpdir, + ) + + output_dir = function_tmpdir / "mf6_output" + result = converter.convert(str(output_dir), verbose=True) + + assert result["grb"].exists(), "GRB file not created" + assert result["hds"].exists(), "MF6 head file not created" + assert result["bud"].exists(), "MF6 budget file not created" + + grb = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + hds_mf6 = HeadFile(str(result["hds"])) + head_mf6 = hds_mf6.get_data(idx=0) + + hds_nwt = HeadFile(str(hds_file)) + head_nwt = hds_nwt.get_data(idx=0) + + np.testing.assert_array_almost_equal( + head_mf6, head_nwt, decimal=5, err_msg="MF6 heads don't match NWT heads" + ) + + bud_mf6 = CellBudgetFile(str(result["bud"])) + texts = bud_mf6.textlist + texts_str = [ + t.decode().strip() if isinstance(t, bytes) else str(t).strip() for t in texts + ] + + assert "FLOW-JA-FACE" in texts_str, "FLOW-JA-FACE not in budget" + assert "DATA-SAT" in texts_str, "DATA-SAT not in budget" + + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", idx=0) + assert flowja is not None, "Could not read FLOW-JA-FACE" + + sat_data = bud_mf6.get_data(text="DATA-SAT", idx=0) + assert sat_data is not None, "Could not read DATA-SAT" + + +@requires_exe("mfnwt") +@pytest.mark.slow +def test_nwt_to_mf6_multilayer_model(function_tmpdir): + from flopy.mf6.utils import MfGrdFile + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowNwt, + ModflowOc, + ModflowRch, + ModflowUpw, + ModflowWel, + ) + from flopy.utils import NwtToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + modelname = "multilayer" + + # 3 layers: top convertible, middle/bottom confined + nlay, nrow, ncol = 3, 10, 10 + delr = delc = 100.0 + + top = np.ones((nrow, ncol)) * 100.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 80.0 + botm[1] = 60.0 + botm[2] = 40.0 + + strt = np.ones((nlay, nrow, ncol)) * 95.0 + + laytyp = np.array([1, 0, 0]) + + hk = np.ones((nlay, nrow, ncol)) * 10.0 + + mf = Modflow( + modelname=modelname, + exe_name="mfnwt", + model_ws=function_tmpdir, + version="mfnwt", + ) + ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + perlen=1, + nstp=1, + steady=True, + ) + ModflowBas(mf, ibound=1, strt=strt) + ModflowUpw(mf, hk=hk, laytyp=laytyp, ipakcb=53) + + # well in center + wel_data = [(1, 5, 5, -1000.0)] # Layer 2, center cell, pumping + ModflowWel(mf, stress_period_data={0: wel_data}) + + # recharge to top layer + ModflowRch(mf, rech=0.001) + + oc = ModflowOc( + mf, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, + ) + oc.reset_budgetunit(budgetunit=53, fname=f"{modelname}.cbc") + + ModflowNwt(mf) + + mf.write_input() + success, _ = mf.run_model(silent=True) + assert success + + hds_file = function_tmpdir / f"{modelname}.hds" + cbc_file = function_tmpdir / f"{modelname}.cbc" + + assert hds_file.exists() + assert cbc_file.exists() + + converter = NwtToMf6Converter( + str(hds_file), + str(cbc_file), + nlay, + nrow, + ncol, + delr=np.ones(ncol) * delr, + delc=np.ones(nrow) * delc, + top=top, + botm=botm, + laytyp=laytyp, + model_ws=function_tmpdir, + ) + + output_dir = function_tmpdir / "mf6_output" + result = converter.convert(str(output_dir), verbose=True) + + # Validate GRB file + grb = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + # Verify ICELLTYPE was set correctly + icelltype = grb._datadict["ICELLTYPE"] + + # Layer 1 should be convertible (1), layers 2-3 confined (0) + # ICELLTYPE is flattened in Fortran order, which interleaves layers + # Pattern: [L0[0,0], L1[0,0], L2[0,0], L0[0,1], L1[0,1], L2[0,1], ...] + ncells = nlay * nrow * ncol + cells_per_layer = nrow * ncol + + # Extract layer data from Fortran-ordered flat array + layer_1_icelltype = icelltype[0::nlay] # Every nlay-th element starting at 0 + layer_2_icelltype = icelltype[1::nlay] # Every nlay-th element starting at 1 + layer_3_icelltype = icelltype[2::nlay] # Every nlay-th element starting at 2 + + assert np.all(layer_1_icelltype == 1), ( + f"Layer 1 should be convertible, got {layer_1_icelltype[:10]}" + ) + assert np.all(layer_2_icelltype == 0), ( + f"Layer 2 should be confined, got {layer_2_icelltype[:10]}" + ) + assert np.all(layer_3_icelltype == 0), ( + f"Layer 3 should be confined, got {layer_3_icelltype[:10]}" + ) + + # Verify TOP array is expanded to all cells + assert len(grb.top) == ncells, f"TOP should have {ncells} values" + + # Extract layer data from Fortran-ordered TOP array + top_layer_1 = grb.top[0::nlay] # Every nlay-th element starting at 0 + top_layer_2 = grb.top[1::nlay] # Every nlay-th element starting at 1 + top_layer_3 = grb.top[2::nlay] # Every nlay-th element starting at 2 + + # Layer 1 top should match model top + np.testing.assert_array_almost_equal( + top_layer_1, + top.flatten(order="F"), + err_msg="Layer 1 TOP should match model top", + ) + + # Layer 2 top should match layer 1 bottom + np.testing.assert_array_almost_equal( + top_layer_2, + botm[0].flatten(order="F"), + err_msg="Layer 2 TOP should match layer 1 bottom", + ) + + # Layer 3 top should match layer 2 bottom + np.testing.assert_array_almost_equal( + top_layer_3, + botm[1].flatten(order="F"), + err_msg="Layer 3 TOP should match layer 2 bottom", + ) + + # Validate heads + hds_nwt = HeadFile(str(hds_file)) + hds_mf6 = HeadFile(str(result["hds"])) + + head_nwt = hds_nwt.get_data(idx=0) + head_mf6 = hds_mf6.get_data(idx=0) + + np.testing.assert_array_almost_equal(head_mf6, head_nwt, decimal=5) + + # Validate budget + bud_mf6 = CellBudgetFile(str(result["bud"])) + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", idx=0) + assert flowja is not None + + # Verify saturation is correct for convertible layer + sat_data = bud_mf6.get_data(text="DATA-SAT", idx=0) + assert sat_data is not None diff --git a/autotest/test_postprocessing.py b/autotest/test_postprocessing.py index d3fbbd2f6f..b2410c2cc7 100644 --- a/autotest/test_postprocessing.py +++ b/autotest/test_postprocessing.py @@ -17,12 +17,17 @@ ModflowTdis, ) from flopy.mf6.utils import get_residuals, get_structured_faceflows +from flopy.mf6.utils.postprocessing import ( + get_structured_connectivity, + get_structured_flowja, +) from flopy.modflow import Modflow, ModflowDis, ModflowLpf, ModflowUpw from flopy.plot import PlotMapView from flopy.utils import get_transmissivities from flopy.utils.gridutil import get_disv_kwargs from flopy.utils.postprocessing import ( get_gradients, + get_saturation, get_specific_discharge, get_water_table, ) @@ -724,3 +729,253 @@ def test_get_transmissivities_fully_saturated(function_tmpdir): assert np.allclose(T_none, T_saturated) assert np.allclose(T_none, T_expected) + + +def test_get_structured_connectivity_simple(): + nlay, nrow, ncol = 2, 2, 2 + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + # Check dimensions + ncells = nlay * nrow * ncol + assert len(ia) == ncells + 1 + assert len(ja) == nja + assert ia[-1] == nja + + # each active cell should have at least 1 connection (diagonal) + for n in range(ncells): + nconn = ia[n + 1] - ia[n] + assert nconn >= 1 + + # first cell (0,0,0) should have 4 connections: diagonal + right + front + lower + nconn_first = ia[1] - ia[0] + first_conns = ja[ia[0] : ia[1]] + assert nconn_first == 4 + assert first_conns[0] == 0 + assert 1 in first_conns + assert 2 in first_conns + assert 4 in first_conns + + +def test_get_structured_connectivity_with_inactive_cells(): + nlay, nrow, ncol = 2, 3, 3 + idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) + idomain[0, 1, 1] = 0 # center cell of top layer inactive + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol, idomain) + + # inactive cell (node 4) should have 0 connections + node_inactive = 0 * nrow * ncol + 1 * ncol + 1 # (k=0, i=1, j=1) + nconn_inactive = ia[node_inactive + 1] - ia[node_inactive] + assert nconn_inactive == 0 + + # neighbors of inactive cell should not connect to it + node_left = 0 * nrow * ncol + 1 * ncol + 0 # (k=0, i=1, j=0) + conns_left = ja[ia[node_left] : ia[node_left + 1]] + assert node_inactive not in conns_left + + +def test_get_structured_connectivity_corner_cells(): + nlay, nrow, ncol = 1, 3, 3 + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + # top-left corner (0,0,0) should have 3 connections: diagonal + right + front + node_corner = 0 + nconn = ia[node_corner + 1] - ia[node_corner] + assert nconn == 3 + + # bottom-right corner (0,2,2) should have 1 connection: diagonal only + node_corner = 0 * nrow * ncol + 2 * ncol + 2 + nconn = ia[node_corner + 1] - ia[node_corner] + assert nconn == 1 + + +def test_get_structured_flowja_to_connections_simple(): + nlay, nrow, ncol = 1, 3, 3 + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + qright = np.ones((nlay, nrow, ncol)) * 1.0 + qfront = np.ones((nlay, nrow, ncol)) * 2.0 + qlower = np.ones((nlay, nrow, ncol)) * 3.0 + flowja = get_structured_flowja( + (qright, qfront, qlower), ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) + + assert len(flowja) == nja, f"flowja length should be {nja}" + + # first cell (0,0,0) connections + node = 0 + conns = ja[ia[node] : ia[node + 1]] + flows = flowja[ia[node] : ia[node + 1]] + for ipos, m in enumerate(conns): + if m == node: + # diagonal + assert flows[ipos] == 0.0 + elif m == 1: + # right + assert flows[ipos] == 1.0 + elif m == 3: + # front + assert flows[ipos] == 2.0 + + +def test_get_structured_flowja_to_connections_multilayer(): + nlay, nrow, ncol = 3, 2, 2 + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + qright = np.ones((nlay, nrow, ncol)) * 10.0 + qfront = np.ones((nlay, nrow, ncol)) * 20.0 + qlower = np.ones((nlay, nrow, ncol)) * 30.0 + flowja = get_structured_flowja( + (qright, qfront, qlower), ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) + + # first cell of top layer should connect to lower layer + node = 0 # (k=0, i=0, j=0) + node_below = nrow * ncol # (k=1, i=0, j=0) + conns = ja[ia[node] : ia[node + 1]] + flows = flowja[ia[node] : ia[node + 1]] + + # check lower connection + lower_idx = np.where(conns == node_below)[0] + if len(lower_idx) > 0: + assert flows[lower_idx[0]] == 30.0 + + +def test_get_saturation_confined_cells(): + nlay, nrow, ncol = 2, 3, 3 + + # all cells confined + icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) + head = np.full((nlay, nrow, ncol), 50.0) + top = np.full((nrow, ncol), 100.0) + botm = np.array( + [ + np.full((nrow, ncol), 75.0), + np.full((nrow, ncol), 25.0), + ] + ) + sat = get_saturation(head, top, botm, icelltype) + + assert np.allclose(sat, 1.0, equal_nan=True) + + +def test_get_saturation_convertible_cells(): + nlay, nrow, ncol = 2, 3, 3 + + # top layer convertible, bottom layer confined + icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) + icelltype[0] = 1 + top = np.full((nrow, ncol), 100.0) + botm = np.array( + [ + np.full((nrow, ncol), 50.0), + np.full((nrow, ncol), 0.0), + ] + ) + head = np.full((nlay, nrow, ncol), 75.0) + sat = get_saturation(head, top, botm, icelltype) + + assert np.allclose(sat[0], 0.5) # unconfined layer + assert np.allclose(sat[1], 1.0) # confined layer + + +def test_get_saturation_fully_saturated(): + nlay, nrow, ncol = 1, 2, 2 + icelltype = np.ones((nlay, nrow, ncol), dtype=np.int32) + top = np.full((nrow, ncol), 100.0) + botm = np.full((nlay, nrow, ncol), 50.0) + head = np.full((nlay, nrow, ncol), 120.0) # above top + sat = get_saturation(head, top, botm, icelltype) + + assert np.allclose(sat, 1.0) # clamped to 1.0 + + +def test_get_saturation_dry_cells(): + nlay, nrow, ncol = 1, 2, 2 + icelltype = np.ones((nlay, nrow, ncol), dtype=np.int32) + top = np.full((nrow, ncol), 100.0) + botm = np.full((nlay, nrow, ncol), 50.0) + head = np.full((nlay, nrow, ncol), 75.0) + head[0, 0, 0] = -999.0 # dry + head[0, 1, 1] = -9999.0 # inactive + sat = get_saturation(head, top, botm, icelltype, hdry=-999.0, hnoflo=-9999.0) + + # dry and inactive cells + assert np.isnan(sat[0, 0, 0]) + assert np.isnan(sat[0, 1, 1]) + + # active cells + assert np.allclose(sat[0, 0, 1], 0.5) + assert np.allclose(sat[0, 1, 0], 0.5) + + +def test_get_saturation_below_bottom(): + nlay, nrow, ncol = 1, 2, 2 + icelltype = np.ones((nlay, nrow, ncol), dtype=np.int32) + top = np.full((nrow, ncol), 100.0) + botm = np.full((nlay, nrow, ncol), 50.0) + head = np.full((nlay, nrow, ncol), 30.0) # below bottom + sat = get_saturation(head, top, botm, icelltype) + + assert np.allclose(sat, 0.0) + + +def test_get_saturation_1d(): + ncells = 10 + icelltype = np.zeros(ncells, dtype=np.int32) + icelltype[0:5] = 1 # first 5 convertible + top = np.linspace(100, 50, ncells) + botm = np.linspace(50, 0, ncells) + head = np.linspace(75, 25, ncells) + sat = get_saturation(head, top, botm, icelltype) + + assert len(sat) == ncells + assert np.allclose(sat[0], 0.5) # convertible cell: (75-50)/(100-50) = 0.5 + assert np.allclose(sat[5:], 1.0) # confined cells: 1.0 + + +def test_get_saturation_invalid_inputs(): + nlay, nrow, ncol = 2, 3, 3 + + head = np.full((nlay, nrow, ncol), 50.0) + top = np.full((nrow, ncol), 100.0) + botm = np.full((nlay, nrow, ncol), 50.0) + icelltype = np.zeros((nlay, nrow - 1, ncol), dtype=np.int32) # Wrong shape + + with pytest.raises(ValueError, match="does not match"): + get_saturation(head, top, botm, icelltype) + + +def test_get_structured_connectivity_big_grid(): + nlay, nrow, ncol = 3, 40, 20 + ncells = nlay * nrow * ncol + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + # Check array sizes + assert len(ia) == ncells + 1 + assert len(ja) == nja + assert ia[-1] == nja + + # Verify all active cells have connections + for n in range(ncells): + nconn = ia[n + 1] - ia[n] + assert nconn >= 1, f"Cell {n} should have at least 1 connection" + + # Interior cell should have 7 connections (diagonal + 6 neighbors) + # But we only store upper triangle, so expect 4 (diagonal + right + front + lower) + node = 1 * nrow * ncol + 10 * ncol + 10 # Middle of layer 1 + nconn = ia[node + 1] - ia[node] + assert nconn == 4, f"Interior cell should have 4 connections, got {nconn}" + + +def test_get_structured_flowja_invalid_inputs(): + nlay, nrow, ncol = 2, 3, 3 + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + qright = np.ones((nlay, nrow, ncol)) + qfront = np.ones((nlay, nrow, ncol)) + qlower = np.ones((nlay, nrow, ncol - 1)) # Wrong shape + + with pytest.raises(ValueError, match="does not match grid shape"): + get_structured_flowja( + (qright, qfront, qlower), ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) diff --git a/flopy/mf6/utils/__init__.py b/flopy/mf6/utils/__init__.py index 45de130dd9..3c527fd78d 100644 --- a/flopy/mf6/utils/__init__.py +++ b/flopy/mf6/utils/__init__.py @@ -3,4 +3,9 @@ from .lakpak_utils import get_lak_connections from .mfsimlistfile import MfSimulationList from .model_splitter import Mf6Splitter -from .postprocessing import get_residuals, get_structured_faceflows +from .postprocessing import ( + get_residuals, + get_structured_connectivity, + get_structured_faceflows, + get_structured_flowja, +) diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index 09c2564543..99f9402dbc 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -931,3 +931,313 @@ def export(self, filename, precision=None, version=1, verbose=False): if verbose: print(f"Successfully wrote {filename}") + + @staticmethod + def write_dis( + filename, + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + ia, + ja, + idomain=None, + icelltype=None, + xorigin=0.0, + yorigin=0.0, + angrot=0.0, + precision="double", + version=1, + verbose=False, + ): + """ + Write a MODFLOW 6 binary grid file (.grb) for a structured (DIS) grid. + + Parameters + ---------- + filename : str or PathLike + Path to output .grb file + nlay : int + Number of layers + nrow : int + Number of rows + ncol : int + Number of columns + delr : array_like + Column spacing, shape (ncol,) + delc : array_like + Row spacing, shape (nrow,) + top : array_like + Top elevation, shape (nrow, ncol) or (ncells,) + botm : array_like + Bottom elevation, shape (nlay, nrow, ncol) or (ncells,) + ia : array_like + CSR row pointers, shape (ncells + 1,), 0-based indexing. + Will be converted to 1-based for the file. + ja : array_like + CSR column indices, shape (nja,), 0-based indexing. + Will be converted to 1-based for the file. + idomain : array_like, optional + Domain array, shape (nlay, nrow, ncol) or (ncells,). + If None, defaults to all active (1). + icelltype : array_like, optional + Cell type array, shape (nlay, nrow, ncol) or (ncells,). + 0 = confined, >0 = convertible. + If None, defaults to all confined (0). + xorigin : float, optional + X-coordinate of grid origin (default 0.0) + yorigin : float, optional + Y-coordinate of grid origin (default 0.0) + angrot : float, optional + Rotation angle in degrees (default 0.0) + precision : str, optional + 'single' or 'double' (default 'double') + version : int, optional + Grid file version (default 1) + verbose : bool, optional + Print progress messages (default False) + + Notes + ----- + The IA and JA arrays should use 0-based indexing (Python convention). + They will be automatically converted to 1-based indexing when written + to the file (Fortran convention). + + Examples + -------- + >>> import numpy as np + >>> from flopy.mf6.utils import MfGrdFile, get_structured_connectivity + >>> nlay, nrow, ncol = 2, 10, 10 + >>> delr = np.ones(ncol) * 100.0 + >>> delc = np.ones(nrow) * 100.0 + >>> top = np.ones((nrow, ncol)) * 10.0 + >>> botm = np.zeros((nlay, nrow, ncol)) + >>> botm[0] = 5.0 + >>> ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + >>> icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) + >>> icelltype[0] = 1 # Top layer convertible + >>> MfGrdFile.write( + ... 'model.dis.grb', + ... nlay, nrow, ncol, + ... delr, delc, top, botm, + ... ia, ja, + ... icelltype=icelltype + ... ) + """ + from pathlib import Path + + import numpy as np + + from ...utils.utils_def import FlopyBinaryData + + # Convert to numpy arrays and handle shapes + delr = np.atleast_1d(delr).astype(np.float64) + delc = np.atleast_1d(delc).astype(np.float64) + ia = np.atleast_1d(ia).astype(np.int32) + ja = np.atleast_1d(ja).astype(np.int32) + + # Handle top - can be (nrow, ncol) or (ncells,) + top = np.asarray(top, dtype=np.float64) + if top.ndim == 2: + if top.shape != (nrow, ncol): + raise ValueError( + f"top shape {top.shape} does not match (nrow, ncol) = ({nrow}, {ncol})" + ) + top = top.flatten(order="F") + elif top.ndim == 1: + # Already flattened + pass + else: + raise ValueError(f"top must be 1D or 2D, got {top.ndim}D") + + # Handle botm - can be (nlay, nrow, ncol) or (ncells,) + botm = np.asarray(botm, dtype=np.float64) + if botm.ndim == 3: + if botm.shape != (nlay, nrow, ncol): + raise ValueError( + f"botm shape {botm.shape} does not match (nlay, nrow, ncol) = ({nlay}, {nrow}, {ncol})" + ) + botm = botm.flatten(order="F") + elif botm.ndim == 1: + # Already flattened + pass + else: + raise ValueError(f"botm must be 1D or 3D, got {botm.ndim}D") + + # Calculate derived values + ncells = nlay * nrow * ncol + nja = len(ja) + + # Set defaults + if idomain is None: + idomain = np.ones(ncells, dtype=np.int32) + else: + idomain = np.atleast_1d(idomain).astype(np.int32) + if idomain.ndim > 1: + idomain = idomain.flatten(order="F") + + if icelltype is None: + icelltype = np.zeros(ncells, dtype=np.int32) + else: + icelltype = np.atleast_1d(icelltype).astype(np.int32) + if icelltype.ndim > 1: + icelltype = icelltype.flatten(order="F") + + # Expand top to all cells if needed + # MF6 expects TOP for every cell (ncells), not just top layer + if len(top) == nrow * ncol: + # Top is model surface only - expand to all layers + # Both TOP and BOTM are in Fortran order (layer-interleaved) + top_all = np.zeros(ncells, dtype=np.float64) + + # For each (row, col) position, set TOP for all layers + # In Fortran order: cells are indexed as (layer, row, col) but + # stored as [L0[0,0], L1[0,0], L2[0,0], L0[0,1], L1[0,1], ...] + for i in range(nrow * ncol): + # Layer 0: use model top + top_all[i * nlay] = top[i] + # Layers 1+: use bottom of layer above + for k in range(1, nlay): + # BOTM is also in Fortran order + # botm[i * nlay + k - 1] = bottom of layer (k-1) at position i + top_all[i * nlay + k] = botm[i * nlay + (k - 1)] + top = top_all + elif len(top) != ncells: + raise ValueError( + f"top length {len(top)} must be nrow*ncol ({nrow * ncol}) or ncells ({ncells})" + ) + + # Validate shapes + if len(delr) != ncol: + raise ValueError(f"delr length {len(delr)} != ncol {ncol}") + if len(delc) != nrow: + raise ValueError(f"delc length {len(delc)} != nrow {nrow}") + if len(botm) != ncells: + raise ValueError(f"botm length {len(botm)} != ncells {ncells}") + if len(ia) != ncells + 1: + raise ValueError(f"ia length {len(ia)} != ncells + 1 ({ncells + 1})") + if len(idomain) != ncells: + raise ValueError(f"idomain length {len(idomain)} != ncells {ncells}") + if len(icelltype) != ncells: + raise ValueError( + f"icelltype length {len(icelltype)} != ncells {ncells}" + ) + + if verbose: + print(f"Writing binary grid file: {filename}") + print(f" Grid type: DIS") + print(f" Dimensions: {nlay} layers, {nrow} rows, {ncol} columns") + print(f" Cells: {ncells}, Connections: {nja}") + print(f" Precision: {precision}") + + # Convert IA/JA to 1-based indexing for file + ia_fortran = ia + 1 + ja_fortran = ja + 1 + + # Build data dictionary + data_dict = { + "NCELLS": ncells, + "NLAY": nlay, + "NROW": nrow, + "NCOL": ncol, + "NJA": nja, + "XORIGIN": xorigin, + "YORIGIN": yorigin, + "ANGROT": angrot, + "DELR": delr, + "DELC": delc, + "TOP": top, + "BOTM": botm, + "IA": ia_fortran, + "JA": ja_fortran, + "IDOMAIN": idomain, + "ICELLTYPE": icelltype, + } + + # Define variable metadata + float_type = "SINGLE" if precision.lower() == "single" else "DOUBLE" + var_list = [ + ("NCELLS", "INTEGER", 0, []), + ("NLAY", "INTEGER", 0, []), + ("NROW", "INTEGER", 0, []), + ("NCOL", "INTEGER", 0, []), + ("NJA", "INTEGER", 0, []), + ("XORIGIN", float_type, 0, []), + ("YORIGIN", float_type, 0, []), + ("ANGROT", float_type, 0, []), + ("DELR", float_type, 1, [ncol]), + ("DELC", float_type, 1, [nrow]), + ("TOP", float_type, 1, [ncells]), + ("BOTM", float_type, 1, [ncells]), + ("IA", "INTEGER", 1, [ncells + 1]), + ("JA", "INTEGER", 1, [nja]), + ("IDOMAIN", "INTEGER", 1, [ncells]), + ("ICELLTYPE", "INTEGER", 1, [ncells]), + ] + + # Create writer with appropriate precision + writer = FlopyBinaryData() + writer.precision = precision + + # Write the file + with open(filename, "wb") as f: + writer.file = f + + # Write header + writer.write_text(f"GRID DIS", nchar=50) + writer.write_text(f"VERSION {version}", nchar=50) + ntxt = len(var_list) + writer.write_text(f"NTXT {ntxt}", nchar=50) + writer.write_text("LENTXT 100", nchar=50) + + # Write variable definitions + for name, dtype_str, ndim, shape in var_list: + if ndim == 0: + # Scalar + defn = f"{name} {dtype_str} NDIM 0" + else: + # Array + shape_str = " ".join(str(s) for s in shape) + defn = f"{name} {dtype_str} NDIM {ndim} {shape_str}" + writer.write_text(defn, nchar=100) + + # Write data + for name, dtype_str, ndim, shape in var_list: + value = data_dict[name] + + if verbose: + if ndim == 0: + print(f" Writing {name} = {value}") + else: + if hasattr(value, "min"): + print( + f" Writing {name}: min = {value.min()}, max = {value.max()}" + ) + else: + print(f" Writing {name}") + + # Write scalar or array data + if ndim == 0: + # Scalar value + if dtype_str == "INTEGER": + writer.write_integer(int(value)) + elif dtype_str in ("DOUBLE", "SINGLE"): + writer.write_real(float(value)) + else: + # Array data + arr = np.asarray(value) + if dtype_str == "INTEGER": + arr = arr.astype(np.int32) + elif dtype_str == "DOUBLE": + arr = arr.astype(np.float64) + elif dtype_str == "SINGLE": + arr = arr.astype(np.float32) + + # Write array (already in correct order from data_dict) + writer.write_record(arr, dtype=arr.dtype) + + if verbose: + print(f"Successfully wrote {filename}") diff --git a/flopy/mf6/utils/postprocessing.py b/flopy/mf6/utils/postprocessing.py index 4f219736d9..47843988fe 100644 --- a/flopy/mf6/utils/postprocessing.py +++ b/flopy/mf6/utils/postprocessing.py @@ -3,6 +3,111 @@ from .binarygrid_util import MfGrdFile +def get_structured_connectivity(nlay, nrow, ncol, idomain=None): + """ + Build IA and JA connectivity arrays for a structured (DIS) grid. + + Parameters + ---------- + nlay : int + Number of layers + nrow : int + Number of rows + ncol : int + Number of columns + idomain : np.ndarray, optional + Domain array indicating active (>0) and inactive (<=0) cells. + Shape: (nlay, nrow, ncol). If None, all cells are active. + + Returns + ------- + ia : np.ndarray + Index array (CSR format), shape (ncells + 1,), dtype int32. + ia[n] is the starting position in ja for cell n's connections. + ia[ncells] is the total number of connections. + ja : np.ndarray + Connection array (CSR format), shape (nja,), dtype int32. + Contains cell numbers for each connection (0-based). + nja : int + Total number of connections + + Notes + ----- + Connectivity order for structured grids (upper triangle only): + 1. Diagonal (self connection) + 2. Right (+1 in j, same k, i) + 3. Front (+1 in i, same k, j) + 4. Lower (+1 in k, same i, j) + + The IA/JA arrays use 0-based indexing (Python convention). + When writing to MF6 binary files, add 1 to convert to Fortran 1-based indexing. + + Examples + -------- + >>> import numpy as np + >>> from flopy.mf6.utils import build_structured_connectivity + >>> nlay, nrow, ncol = 2, 3, 3 + >>> ia, ja, nja = build_structured_connectivity(nlay, nrow, ncol) + >>> print(f"Total cells: {nlay * nrow * ncol}, connections: {nja}") + Total cells: 18, connections: 42 + """ + ncells = nlay * nrow * ncol + + # Default to all active cells if idomain not provided + if idomain is None: + idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) + else: + idomain = np.asarray(idomain, dtype=np.int32) + if idomain.shape != (nlay, nrow, ncol): + raise ValueError( + f"idomain shape {idomain.shape} does not match grid shape " + f"({nlay}, {nrow}, {ncol})" + ) + + ia = np.zeros(ncells + 1, dtype=np.int32) + ja_list = [] + nja = 0 + + for k in range(nlay): + for i in range(nrow): + for j in range(ncol): + node = k * nrow * ncol + i * ncol + j + + # Skip inactive cells - they still get an entry in IA + if idomain[k, i, j] <= 0: + ia[node + 1] = nja + continue + + # Add diagonal (self connection) + ja_list.append(node) + nja += 1 + + # Add connections to neighbors (upper triangle only) + # Right neighbor (j+1) + if j + 1 < ncol and idomain[k, i, j + 1] > 0: + m = k * nrow * ncol + i * ncol + (j + 1) + ja_list.append(m) + nja += 1 + + # Front neighbor (i+1) + if i + 1 < nrow and idomain[k, i + 1, j] > 0: + m = k * nrow * ncol + (i + 1) * ncol + j + ja_list.append(m) + nja += 1 + + # Lower neighbor (k+1) + if k + 1 < nlay and idomain[k + 1, i, j] > 0: + m = (k + 1) * nrow * ncol + i * ncol + j + ja_list.append(m) + nja += 1 + + ia[node + 1] = nja + + ja = np.array(ja_list, dtype=np.int32) + + return ia, ja, nja + + def get_structured_faceflows( flowja, grb_file=None, @@ -52,19 +157,12 @@ def get_structured_faceflows( grb = MfGrdFile(grb_file, verbose=verbose) if grb.grid_type != "DIS": raise ValueError( - "get_structured_faceflows method " - "is only for structured DIS grids" + "get_structured_faceflows method is only for structured DIS grids" ) ia, ja = grb.ia, grb.ja nlay, nrow, ncol = grb.nlay, grb.nrow, grb.ncol else: - if ( - ia is None - or ja is None - or nlay is None - or nrow is None - or ncol is None - ): + if ia is None or ja is None or nlay is None or nrow is None or ncol is None: raise ValueError( "ia, ja, nlay, nrow, and ncol must be" "specified if a MODFLOW 6 binary grid" @@ -76,7 +174,7 @@ def get_structured_faceflows( flowja = flowja.flatten() # evaluate size of flowja relative to ja - __check_flowja_size(flowja, ja) + _check_flowja_size(flowja, ja) # create empty flat face flow arrays shape = (nlay, nrow, ncol) @@ -134,7 +232,8 @@ def get_face(m, n, nlay, nrow, ncol): # fill right, front and lower face flows # (below main diagonal) flows = [frf, fff, flf] - for n in range(grb.nodes): + nodes = nlay * nrow * ncol + for n in range(nodes): for i in range(ia[n] + 1, ia[n + 1]): m = ja[i] if m <= n: @@ -146,9 +245,147 @@ def get_face(m, n, nlay, nrow, ncol): return frf.reshape(shape), fff.reshape(shape), flf.reshape(shape) -def get_residuals( - flowja, grb_file=None, ia=None, ja=None, shape=None, verbose=False +def get_structured_flowja( + faceflows, + grb_file=None, + ia=None, + ja=None, + nlay=None, + nrow=None, + ncol=None, + idomain=None, + verbose=False, ): + """ + Get connection flows (flowja) from face flows for a structured grid. + + This is the inverse of get_structured_faceflows(). Converts MODFLOW-2005/NWT + style face flows (flow right face, flow front face, flow lower face) to + MODFLOW 6 style connection flows for the FLOW-JA-FACE budget term. + + Parameters + ---------- + faceflows : tuple of ndarray + Tuple of (frf, fff, flf) where: + - frf : flow right face, shape (nlay, nrow, ncol) + - fff : flow front face, shape (nlay, nrow, ncol) + - flf : flow lower face, shape (nlay, nrow, ncol) + grb_file : str, optional + MODFLOW 6 binary grid file path + ia : list or ndarray, optional + CRS row pointers. Only required if grb_file is not provided. + ja : list or ndarray, optional + CRS column pointers. Only required if grb_file is not provided. + nlay : int, optional + Number of layers. Only required if grb_file is not provided. + nrow : int, optional + Number of rows. Only required if grb_file is not provided. + ncol : int, optional + Number of columns. Only required if grb_file is not provided. + idomain : ndarray, optional + Domain array, shape (nlay, nrow, ncol) + verbose : bool, optional + Write information to standard output (default False) + + Returns + ------- + flowja : ndarray + Connection flows, size (nja,) + + See Also + -------- + get_structured_faceflows : Inverse operation (flowja to face flows) + + Examples + -------- + >>> import numpy as np + >>> from flopy.mf6.utils import get_structured_flowja + >>> nlay, nrow, ncol = 1, 3, 3 + >>> frf = np.ones((nlay, nrow, ncol)) * 1.0 + >>> fff = np.ones((nlay, nrow, ncol)) * 2.0 + >>> flf = np.ones((nlay, nrow, ncol)) * 3.0 + >>> flowja = get_structured_flowja((frf, fff, flf), + ... nlay=nlay, nrow=nrow, ncol=ncol) + """ + # Unpack face flows + qright, qfront, qlower = faceflows + + # Get grid information + if grb_file is not None: + grb = MfGrdFile(grb_file, verbose=verbose) + if grb.grid_type != "DIS": + raise ValueError( + "get_structured_flowja method is only for structured DIS grids" + ) + ia, ja = grb.ia, grb.ja + nlay, nrow, ncol = grb.nlay, grb.nrow, grb.ncol + else: + if ia is None or ja is None or nlay is None or nrow is None or ncol is None: + raise ValueError( + "ia, ja, nlay, nrow, and ncol must be specified if grb_file is not provided" + ) + + # Validate input shapes + expected_shape = (nlay, nrow, ncol) + for name, arr in [("qright", qright), ("qfront", qfront), ("qlower", qlower)]: + arr = np.asarray(arr) + if arr.shape != expected_shape: + raise ValueError( + f"{name} shape {arr.shape} does not match grid shape {expected_shape}" + ) + + # Convert to arrays + qright = np.asarray(qright, dtype=np.float64) + qfront = np.asarray(qfront, dtype=np.float64) + qlower = np.asarray(qlower, dtype=np.float64) + + # Default to all active if idomain not provided + if idomain is None: + idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) + else: + idomain = np.asarray(idomain, dtype=np.int32) + + ncells = nlay * nrow * ncol + nja = len(ja) + flowja = np.zeros(nja, dtype=np.float64) + + for n in range(ncells): + # Skip inactive cells + k, i, j = np.unravel_index(n, (nlay, nrow, ncol)) + if idomain[k, i, j] <= 0: + continue + + # Get connections for this cell + istart = ia[n] + iend = ia[n + 1] + + for ipos in range(istart, iend): + m = ja[ipos] + + # Diagonal - no self flow + if m == n: + flowja[ipos] = 0.0 + continue + + # Determine connection type by comparing node numbers + km, im, jm = np.unravel_index(m, (nlay, nrow, ncol)) + + # Right connection (j increases by 1) + if km == k and im == i and jm == j + 1: + flowja[ipos] = qright[k, i, j] + + # Front connection (i increases by 1) + elif km == k and im == i + 1 and jm == j: + flowja[ipos] = qfront[k, i, j] + + # Lower connection (k increases by 1) + elif km == k + 1 and im == i and jm == j: + flowja[ipos] = qlower[k, i, j] + + return flowja + + +def get_residuals(flowja, grb_file=None, ia=None, ja=None, shape=None, verbose=False): """ Get the residual from the MODFLOW 6 flowja flows. The residual is stored in the diagonal position of the flowja vector. @@ -191,7 +428,7 @@ def get_residuals( flowja = flowja.flatten() # evaluate size of flowja relative to ja - __check_flowja_size(flowja, ja) + _check_flowja_size(flowja, ja) # create residual nodes = grb.nodes @@ -211,12 +448,9 @@ def get_residuals( return residual -# internal -def __check_flowja_size(flowja, ja): +def _check_flowja_size(flowja, ja): """ Check the shape of flowja relative to ja. """ if flowja.shape != ja.shape: - raise ValueError( - f"size of flowja ({flowja.shape}) not equal to {ja.shape}" - ) + raise ValueError(f"size of flowja ({flowja.shape}) not equal to {ja.shape}") diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index bd15a59b86..5d628fc860 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -40,6 +40,7 @@ from .mfreadnam import parsenamefile from .modpathfile import EndpointFile, PathlineFile, TimeseriesFile from .mtlistfile import MtListBudget +from .nwt_to_mf6 import NwtToMf6Converter, get_icelltype_from_laytyp from .observationfile import HydmodObs, Mf6Obs, SwrObs from .optionblock import OptionBlock from .postprocessing import get_specific_discharge, get_transmissivities diff --git a/flopy/utils/nwt_to_mf6.py b/flopy/utils/nwt_to_mf6.py new file mode 100644 index 0000000000..8ac5c0149b --- /dev/null +++ b/flopy/utils/nwt_to_mf6.py @@ -0,0 +1,521 @@ +""" +Utilities for converting MODFLOW-NWT binary outputs to MODFLOW 6 format. +""" + +import os +from pathlib import Path + +import numpy as np + + +def get_icelltype_from_laytyp(laytyp): + """ + Convert NWT LAYTYP values to MF6 ICELLTYPE values. + + Parameters + ---------- + laytyp : array_like + Layer type array from NWT LPF or UPW package. + - 0: Confined + - >0: Convertible (unconfined/water table) + Shape can be (nlay,) or (nlay, nrow, ncol) + + Returns + ------- + icelltype : ndarray + Cell type array for MF6 (same shape as input): + - 0: Confined + - 1: Convertible + + Notes + ----- + In MODFLOW-NWT, LAYTYP indicates layer properties: + - LAYTYP = 0: Confined, transmissivity constant + - LAYTYP > 0: Convertible, transmissivity varies with saturation + + In MODFLOW 6, ICELLTYPE similarly indicates cell properties: + - ICELLTYPE = 0: Confined + - ICELLTYPE > 0: Convertible + + For conversion, we use a simple mapping: + - LAYTYP = 0 → ICELLTYPE = 0 + - LAYTYP > 0 → ICELLTYPE = 1 + + Examples + -------- + >>> import numpy as np + >>> from flopy.utils.nwt_to_mf6 import get_icelltype_from_laytyp + >>> laytyp = np.array([1, 0, 0]) # Top layer convertible + >>> icelltype = get_icelltype_from_laytyp(laytyp) + >>> print(icelltype) + [1 0 0] + """ + laytyp = np.atleast_1d(laytyp).astype(np.int32) + icelltype = np.where(laytyp > 0, 1, 0).astype(np.int32) + return icelltype + + +class NwtToMf6Converter: + """ + Convert MODFLOW-NWT binary outputs to MODFLOW 6 binary format. + + This class reads NWT head and budget files, applies necessary + transformations, and writes MF6-compatible binary files that can + be consumed by PRT or other MF6 post-processors. + + Parameters + ---------- + hds_file : str + Path to NWT head file (.hds) + cbc_file : str + Path to NWT cell budget file (.cbc) + nlay : int + Number of layers + nrow : int + Number of rows + ncol : int + Number of columns + delr : array_like + Column spacing, shape (ncol,) + delc : array_like + Row spacing, shape (nrow,) + top : array_like + Top elevation, shape (nrow, ncol) + botm : array_like + Bottom elevation, shape (nlay, nrow, ncol) + laytyp : array_like + Layer type from LPF/UPW, shape (nlay,). + 0 = confined, >0 = convertible + idomain : array_like, optional + Domain array, shape (nlay, nrow, ncol). + >0 = active, 0 = inactive, <0 = vertical pass-through + If None, all cells are active. + hdry : float, optional + Head value for dry cells (default -999.0) + hnoflo : float, optional + Head value for inactive cells (default -9999.0) + model_ws : str or PathLike, optional + Model workspace for input files (default current directory) + + Examples + -------- + >>> import numpy as np + >>> from flopy.utils import NwtToMf6Converter + >>> # Set up grid parameters + >>> nlay, nrow, ncol = 3, 10, 10 + >>> delr = np.ones(ncol) * 100.0 + >>> delc = np.ones(nrow) * 100.0 + >>> top = np.ones((nrow, ncol)) * 10.0 + >>> botm = np.zeros((nlay, nrow, ncol)) + >>> botm[0] = 5.0 + >>> botm[1] = 0.0 + >>> botm[2] = -5.0 + >>> laytyp = np.array([1, 0, 0]) # Top layer convertible + >>> + >>> # Create converter + >>> converter = NwtToMf6Converter( + ... 'model.hds', + ... 'model.cbc', + ... nlay, nrow, ncol, + ... delr, delc, top, botm, + ... laytyp + ... ) + >>> + >>> # Convert files + >>> converter.convert('mf6_output') + """ + + def __init__( + self, + hds_file, + cbc_file, + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + laytyp, + idomain=None, + hdry=-999.0, + hnoflo=-9999.0, + model_ws=".", + ): + from ..mf6.utils import get_structured_connectivity + from .binaryfile import CellBudgetFile, HeadFile + + self.hds_file = Path(model_ws) / hds_file + self.cbc_file = Path(model_ws) / cbc_file + self.nlay = nlay + self.nrow = nrow + self.ncol = ncol + self.ncells = nlay * nrow * ncol + + # Store grid geometry + self.delr = np.atleast_1d(delr).astype(np.float64) + self.delc = np.atleast_1d(delc).astype(np.float64) + self.top = np.atleast_2d(top).astype(np.float64) + self.botm = np.atleast_3d(botm).astype(np.float64) + + # Convert LAYTYP to ICELLTYPE + self.laytyp = np.atleast_1d(laytyp).astype(np.int32) + self.icelltype = get_icelltype_from_laytyp(laytyp) + + # Expand ICELLTYPE to 3D if needed + if self.icelltype.ndim == 1: + # Expand (nlay,) to (nlay, nrow, ncol) + # Use broadcasting to properly replicate each layer's value + self.icelltype_3d = np.broadcast_to( + self.icelltype[:, np.newaxis, np.newaxis], + (nlay, nrow, ncol), + ).copy() # Copy to make it writable + else: + self.icelltype_3d = self.icelltype + + # Set IDOMAIN + if idomain is None: + self.idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) + else: + self.idomain = np.atleast_3d(idomain).astype(np.int32) + + self.hdry = hdry + self.hnoflo = hnoflo + + # Build connectivity arrays + self.ia, self.ja, self.nja = get_structured_connectivity( + nlay, nrow, ncol, self.idomain + ) + + # Open binary files + self.hds_obj = HeadFile(str(self.hds_file)) + self.cbc_obj = CellBudgetFile(str(self.cbc_file)) + + # Get time information + self.times = self.hds_obj.get_times() + self.kstpkper = self.hds_obj.get_kstpkper() + + def convert( + self, + output_dir, + grb_name="gwf.grb", + hds_name="gwf.hds", + bud_name="gwf.bud", + precision="double", + verbose=False, + ): + """ + Convert NWT binary outputs to MF6 format. + + Parameters + ---------- + output_dir : str or PathLike + Directory for output files (will be created if doesn't exist) + grb_name : str, optional + Grid file name (default 'gwf.grb') + hds_name : str, optional + Head file name (default 'gwf.hds') + bud_name : str, optional + Budget file name (default 'gwf.bud') + precision : str, optional + 'single' or 'double' for binary files (default 'double') + verbose : bool, optional + Print progress messages (default False) + + Returns + ------- + dict + Paths to created files: + - 'grb': Path to grid file + - 'hds': Path to head file + - 'bud': Path to budget file + """ + from ..mf6.utils import MfGrdFile + from .binaryfile import CellBudgetFile, HeadFile + + # Create output directory + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + grb_path = output_dir / grb_name + hds_path = output_dir / hds_name + bud_path = output_dir / bud_name + + if verbose: + print("\nConverting NWT outputs to MF6 format") + print(f" Input head file: {self.hds_file}") + print(f" Input budget file: {self.cbc_file}") + print(f" Output directory: {output_dir}") + print(f" Grid dimensions: {self.nlay}L x {self.nrow}R x {self.ncol}C") + print(f" Time steps: {len(self.times)}") + + # Write GRB file + if verbose: + print(f"\nWriting grid file: {grb_path}") + MfGrdFile.write_dis( + str(grb_path), + self.nlay, + self.nrow, + self.ncol, + self.delr, + self.delc, + self.top, + self.botm, + self.ia, + self.ja, + idomain=self.idomain, + icelltype=self.icelltype_3d, + precision=precision, + verbose=verbose, + ) + + # Write HDS file + if verbose: + print(f"\nWriting head file: {hds_path}") + self._write_heads(hds_path, precision, verbose) + + # Write BUD file + if verbose: + print(f"\nWriting budget file: {bud_path}") + self._write_budget(bud_path, precision, verbose) + + if verbose: + print("\nConversion complete!") + + return {"grb": grb_path, "hds": hds_path, "bud": bud_path} + + def _write_heads(self, filename, precision, verbose): + """Write head file with all time steps.""" + from .binaryfile import HeadFile + + # Read all heads + head_dict = {} + totim_dict = {} + pertim_dict = {} + + for idx, kstpkper in enumerate(self.kstpkper): + head = self.hds_obj.get_data(kstpkper=kstpkper) + totim = self.times[idx] + # Assume pertim = totim for now (could be refined) + pertim = totim + + head_dict[kstpkper] = head + totim_dict[kstpkper] = totim + pertim_dict[kstpkper] = pertim + + if verbose: + print( + f" Read head for time step {kstpkper}: " + f"min={head.min():.2f}, max={head.max():.2f}" + ) + + # Write using HeadFile.write() + HeadFile.write( + str(filename), + head_dict, + nlay=self.nlay, + nrow=self.nrow, + ncol=self.ncol, + precision=precision, + totim=totim_dict, + pertim=pertim_dict, + verbose=verbose, + ) + + def _write_budget(self, filename, precision, verbose): + """Write budget file with FLOW-JA-FACE, DATA-SPDIS, and DATA-SAT.""" + from ..mf6.utils import get_structured_flowja + from .binaryfile import CellBudgetFile + from .postprocessing import get_saturation + + if verbose: + print(f" Processing {len(self.kstpkper)} time steps...") + + # We'll write three terms per time step: + # 1. FLOW-JA-FACE (imeth=1, array) + # 2. DATA-SPDIS (imeth=6, list with qx, qy, qz) + # 3. DATA-SAT (imeth=6, list) + + # Build list of records + records = [] + + for idx, kstpkper in enumerate(self.kstpkper): + kstp, kper = kstpkper + totim = self.times[idx] + pertim = totim # Simplified + delt = 1.0 # Simplified - could calculate from times + + if verbose: + print(f" Processing time step {kstpkper}...") + + # Get head + head = self.hds_obj.get_data(kstpkper=kstpkper) + + # Get face flows + if verbose: + print( + f" Available budget terms: " + f"{self.cbc_obj.get_unique_record_names()}" + ) + + # Check which face flows are available + # For 1D/2D models, not all face flows may exist + available_terms = [ + t.decode().strip() for t in self.cbc_obj.get_unique_record_names() + ] + + try: + # FLOW RIGHT FACE (required for X-direction flow) + if "FLOW RIGHT FACE" in available_terms: + frf_data = self.cbc_obj.get_data( + text="FLOW RIGHT FACE", kstpkper=kstpkper + ) + if verbose: + print( + f" FLOW RIGHT FACE: {type(frf_data)}, " + f"len={len(frf_data) if frf_data else 0}" + ) + frf = frf_data[0] if frf_data and len(frf_data) > 0 else None + else: + frf = None + + # FLOW FRONT FACE (required for Y-direction flow) + if "FLOW FRONT FACE" in available_terms: + fff_data = self.cbc_obj.get_data( + text="FLOW FRONT FACE", kstpkper=kstpkper + ) + if verbose: + print( + f" FLOW FRONT FACE: {type(fff_data)}, " + f"len={len(fff_data) if fff_data else 0}" + ) + fff = fff_data[0] if fff_data and len(fff_data) > 0 else None + else: + fff = None + + # FLOW LOWER FACE (required for Z-direction flow) + if "FLOW LOWER FACE" in available_terms: + flf_data = self.cbc_obj.get_data( + text="FLOW LOWER FACE", kstpkper=kstpkper + ) + if verbose: + print( + f" FLOW LOWER FACE: {type(flf_data)}, " + f"len={len(flf_data) if flf_data else 0}" + ) + flf = flf_data[0] if flf_data and len(flf_data) > 0 else None + else: + flf = None + + # Validate at least one face flow exists + if frf is None and fff is None and flf is None: + raise ValueError("No face flows found in budget file") + + # Create zero arrays for missing face flows + # For 1D/2D models, not all directions have flow + shape_3d = (self.nlay, self.nrow, self.ncol) + if frf is None: + frf = np.zeros(shape_3d, dtype=np.float64) + if fff is None: + fff = np.zeros(shape_3d, dtype=np.float64) + if flf is None: + flf = np.zeros(shape_3d, dtype=np.float64) + + except Exception as e: + if verbose: + print(f" Warning: Could not read face flows: {e}") + print(f" Skipping time step {kstpkper}") + continue + + # 1. Convert to FLOW-JA-FACE + flowja = get_structured_flowja( + (frf, fff, flf), + ia=self.ia, + ja=self.ja, + nlay=self.nlay, + nrow=self.nrow, + ncol=self.ncol, + ) + + records.append( + { + "data": flowja, + "kstp": kstp, + "kper": kper, + "totim": totim, + "pertim": pertim, + "delt": delt, + "text": "FLOW-JA-FACE", + "imeth": 1, + } + ) + + # 2. DATA-SPDIS (specific discharge) - SKIPPED for now + # TODO: Requires refactoring get_specific_discharge() to work + # without model object or implementing a minimal wrapper. PRT can + # calculate specific discharge from FLOW-JA-FACE if needed. + # See PHASE3_REMAINING_ISSUES.md for details. + + # 3. Calculate saturation + sat = get_saturation( + head, self.top, self.botm, self.icelltype_3d, self.hdry, self.hnoflo + ) + + # Build list data for DATA-SAT + sat_flat = sat.flatten(order="F") + active_sat = ~np.isnan(sat_flat) + nlist_sat = np.sum(active_sat) + + if nlist_sat > 0: + nodes_sat = np.arange(self.ncells)[active_sat] + 1 # 1-based + + # Create structured array for imeth=6 + dtype = np.dtype( + [ + ("node", np.int32), + ("node2", np.int32), + ("sat", np.float64), + ] + ) + sat_data = np.zeros(nlist_sat, dtype=dtype) + sat_data["node"] = nodes_sat + sat_data["node2"] = nodes_sat + sat_data["sat"] = sat_flat[active_sat] + + records.append( + { + "data": sat_data, + "kstp": kstp, + "kper": kper, + "totim": totim, + "pertim": pertim, + "delt": delt, + "text": "DATA-SAT", + "imeth": 6, + "ndat": 1, + } + ) + + # Write all records + if verbose: + print(f" Writing {len(records)} budget records...") + + CellBudgetFile.write( + str(filename), + records, + precision=precision, + nlay=self.nlay, + nrow=self.nrow, + ncol=self.ncol, + verbose=verbose, + ) + + def __repr__(self): + return ( + f"NwtToMf6Converter(\n" + f" hds_file={self.hds_file},\n" + f" cbc_file={self.cbc_file},\n" + f" grid={self.nlay}x{self.nrow}x{self.ncol},\n" + f" time_steps={len(self.times)}\n" + f")" + ) diff --git a/flopy/utils/postprocessing.py b/flopy/utils/postprocessing.py index ec94166267..d393ef034c 100644 --- a/flopy/utils/postprocessing.py +++ b/flopy/utils/postprocessing.py @@ -939,3 +939,200 @@ def get_specific_discharge( qz[noflo_or_dry] = np.nan return qx, qy, qz + + +def get_saturation(head, top, botm, icelltype, hdry=-999.0, hnoflo=-9999.0): + """ + Calculate cell saturation from head values. + + Computes the fraction of each cell that is saturated based on head, + cell top/bottom elevations, and cell type. + + Parameters + ---------- + head : np.ndarray + Head values, shape (nlay, nrow, ncol) or (ncells,) + top : np.ndarray + Top elevation of cells, shape (nrow, ncol) for top layer or + (nlay, nrow, ncol) for all layers, or (ncells,) + botm : np.ndarray + Bottom elevation of cells, shape (nlay, nrow, ncol) or (ncells,) + icelltype : np.ndarray + Cell type indicator: + - 0: confined (always fully saturated) + - >0: convertible/unconfined (saturation varies with head) + Shape: (nlay, nrow, ncol) or (ncells,) + hdry : float, optional + Head value indicating dry cell (default -999.0) + hnoflo : float, optional + Head value indicating inactive cell (default -9999.0) + + Returns + ------- + saturation : np.ndarray + Cell saturation values (0.0 to 1.0), same shape as head. + - 1.0 for fully saturated cells (confined or head >= top) + - 0.0 to 1.0 for partially saturated cells + - NaN for inactive or dry cells + + Notes + ----- + Saturation calculation: + + For confined cells (icelltype == 0): + sat = 1.0 + + For convertible cells (icelltype > 0): + sat = (head - botm) / (top - botm) + Clamped to [0.0, 1.0] + + For dry or inactive cells: + sat = NaN + + The top elevation for layer k is: + - Layer 0: top array (model top) + - Layer k>0: botm[k-1] (bottom of layer above) + + Examples + -------- + >>> import numpy as np + >>> from flopy.utils.postprocessing import calculate_saturation + >>> nlay, nrow, ncol = 3, 10, 10 + >>> head = np.full((nlay, nrow, ncol), 50.0) + >>> top = np.full((nrow, ncol), 100.0) + >>> botm = np.array([ + ... np.full((nrow, ncol), 75.0), + ... np.full((nrow, ncol), 50.0), + ... np.full((nrow, ncol), 25.0) + ... ]) + >>> icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) + >>> icelltype[0] = 1 # Top layer convertible + >>> sat = calculate_saturation(head, top, botm, icelltype) + >>> print(f"Layer 0 saturation: {sat[0,0,0]:.2f}") # Partially saturated + Layer 0 saturation: 0.00 + >>> print(f"Layer 1 saturation: {sat[1,0,0]:.2f}") # Fully saturated (confined) + Layer 1 saturation: 1.00 + """ + # Convert to arrays + head = np.asarray(head, dtype=np.float64) + top = np.asarray(top, dtype=np.float64) + botm = np.asarray(botm, dtype=np.float64) + icelltype = np.asarray(icelltype, dtype=np.int32) + + # Determine if arrays are 3D or 1D + is_3d = head.ndim == 3 + + if is_3d: + nlay, nrow, ncol = head.shape + ncells = nlay * nrow * ncol + + # Ensure top is 2D for structured grids + if top.ndim == 3: + # If top is 3D, use first layer + top2d = top[0] + elif top.ndim == 2: + top2d = top + else: + raise ValueError( + f"top must be 2D (nrow, ncol) or 3D (nlay, nrow, ncol), " + f"got shape {top.shape}" + ) + + # Ensure botm is 3D + if botm.ndim == 3: + botm3d = botm + else: + raise ValueError( + f"botm must be 3D (nlay, nrow, ncol), got shape {botm.shape}" + ) + + # Ensure icelltype matches head shape + if icelltype.shape != head.shape: + raise ValueError( + f"icelltype shape {icelltype.shape} does not match " + f"head shape {head.shape}" + ) + + # Initialize saturation + sat = np.ones_like(head, dtype=np.float64) + + # Process each cell + for k in range(nlay): + for i in range(nrow): + for j in range(ncol): + h = head[k, i, j] + + # Check for inactive or dry cells + if h <= hnoflo or h <= hdry: + sat[k, i, j] = np.nan + continue + + # Confined cells are always fully saturated + if icelltype[k, i, j] == 0: + sat[k, i, j] = 1.0 + continue + + # Convertible cell - calculate saturation + if k == 0: + top_elev = top2d[i, j] + else: + top_elev = botm3d[k - 1, i, j] + + bot_elev = botm3d[k, i, j] + thickness = top_elev - bot_elev + + if thickness <= 0: + # Invalid cell geometry + sat[k, i, j] = np.nan + continue + + # Calculate saturated thickness + sat_thickness = max(0.0, min(h - bot_elev, thickness)) + sat[k, i, j] = sat_thickness / thickness + + else: + # 1D arrays (unstructured or flattened) + ncells = len(head) + + # All arrays must be 1D + if top.ndim != 1 or botm.ndim != 1 or icelltype.ndim != 1: + raise ValueError("For 1D head, top, botm, and icelltype must all be 1D") + + if len(top) != ncells or len(botm) != ncells or len(icelltype) != ncells: + raise ValueError( + f"All arrays must have same length: head={ncells}, " + f"top={len(top)}, botm={len(botm)}, icelltype={len(icelltype)}" + ) + + # Initialize saturation + sat = np.ones(ncells, dtype=np.float64) + + # Process each cell + for n in range(ncells): + h = head[n] + + # Check for inactive or dry cells + if h <= hnoflo or h <= hdry: + sat[n] = np.nan + continue + + # Confined cells are always fully saturated + if icelltype[n] == 0: + sat[n] = 1.0 + continue + + # Convertible cell - calculate saturation + top_elev = top[n] + bot_elev = botm[n] + thickness = top_elev - bot_elev + + if thickness <= 0: + # Invalid cell geometry + sat[n] = np.nan + continue + + # Calculate saturated thickness + sat_thickness = max(0.0, min(h - bot_elev, thickness)) + sat[n] = sat_thickness / thickness + + return sat From 39eeb101131bcc19008a65aeef1e598a075f7312 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 15 Apr 2026 11:53:13 -0700 Subject: [PATCH 02/10] negative laytyp means convertible, bidirectional connectivity, pertim/delt fix --- autotest/test_nwt_to_mf6.py | 11 ++- autotest/test_postprocessing.py | 80 +++++++++++++----- flopy/mf6/utils/postprocessing.py | 133 +++++++++++++++++------------- flopy/utils/nwt_to_mf6.py | 65 +++++++++++---- 4 files changed, 194 insertions(+), 95 deletions(-) diff --git a/autotest/test_nwt_to_mf6.py b/autotest/test_nwt_to_mf6.py index 7376f1f2ee..650201d17e 100644 --- a/autotest/test_nwt_to_mf6.py +++ b/autotest/test_nwt_to_mf6.py @@ -21,11 +21,16 @@ def test_get_icelltype_from_laytyp(): icelltype = get_icelltype_from_laytyp(laytyp) assert icelltype == 1 - # various LAYTYP values + # negative laytyp is convertible (wetting disabled, not confined) + laytyp = -1 + icelltype = get_icelltype_from_laytyp(laytyp) + assert icelltype == 1 + + # various LAYTYP values — any non-zero means convertible laytyp = np.array([0, 1, 2, 3, -1]) icelltype = get_icelltype_from_laytyp(laytyp) - # All > 0 map to 1, 0 stays 0, <0 stays 0 - assert np.array_equal(icelltype, [0, 1, 1, 1, 0]) + # 0 stays 0, all others (positive or negative) map to 1 + assert np.array_equal(icelltype, [0, 1, 1, 1, 1]) def test_mfgrdfile_write_roundtrip(tmp_path): diff --git a/autotest/test_postprocessing.py b/autotest/test_postprocessing.py index b2410c2cc7..41cb8c5f5a 100644 --- a/autotest/test_postprocessing.py +++ b/autotest/test_postprocessing.py @@ -777,15 +777,22 @@ def test_get_structured_connectivity_corner_cells(): nlay, nrow, ncol = 1, 3, 3 ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) - # top-left corner (0,0,0) should have 3 connections: diagonal + right + front + # top-left corner (0,0,0): no left or back neighbor at boundary + # connections: diagonal + right + front = 3 node_corner = 0 nconn = ia[node_corner + 1] - ia[node_corner] assert nconn == 3 - # bottom-right corner (0,2,2) should have 1 connection: diagonal only + # bottom-right corner (0,2,2): no right or front neighbor at boundary + # connections: diagonal + left + back = 3 (full bidirectional connectivity) node_corner = 0 * nrow * ncol + 2 * ncol + 2 - nconn = ia[node_corner + 1] - ia[node_corner] - assert nconn == 1 + conns = ja[ia[node_corner] : ia[node_corner + 1]] + nconn = len(conns) + assert nconn == 3 + # left neighbor (0,2,1) = node 7, back neighbor (0,1,2) = node 5 + assert node_corner in conns # diagonal + assert 7 in conns # left + assert 5 in conns # back def test_get_structured_flowja_to_connections_simple(): @@ -800,20 +807,36 @@ def test_get_structured_flowja_to_connections_simple(): assert len(flowja) == nja, f"flowja length should be {nja}" - # first cell (0,0,0) connections + # Corner cell (0,0,0): only right and front outflow, no lower (1 layer) + # MODFLOW 6 sign convention: outflow from n is negative, inflow is positive. node = 0 conns = ja[ia[node] : ia[node + 1]] flows = flowja[ia[node] : ia[node + 1]] - for ipos, m in enumerate(conns): - if m == node: - # diagonal - assert flows[ipos] == 0.0 - elif m == 1: - # right - assert flows[ipos] == 1.0 - elif m == 3: - # front - assert flows[ipos] == 2.0 + conn_map = dict(zip(conns, flows)) + assert conn_map[node] == 0.0 # diagonal + assert conn_map[1] == -1.0 # right outflow: -qright[0,0,0] + assert conn_map[3] == -2.0 # front outflow: -qfront[0,0,0] + + # Interior cell (0,1,1): left/back/right/front all present + node = 1 * ncol + 1 # (k=0, i=1, j=1) + conns = ja[ia[node] : ia[node + 1]] + flows = flowja[ia[node] : ia[node + 1]] + conn_map = dict(zip(conns, flows)) + assert conn_map[node] == 0.0 # diagonal + assert conn_map[node - 1] == qright[0, 1, 0] # left inflow: +qright[0,1,0] + assert conn_map[node + 1] == -qright[0, 1, 1] # right outflow: -qright[0,1,1] + assert conn_map[node - ncol] == qfront[0, 0, 1] # back inflow: +qfront[0,0,1] + assert conn_map[node + ncol] == -qfront[0, 1, 1] # front outflow: -qfront[0,1,1] + + # Verify roundtrip: flowja → get_structured_faceflows → should recover originals. + # Boundary cells with no neighbor in that direction have no stored connection, + # so the roundtrip can only recover values for interior faces. + from flopy.mf6.utils import get_structured_faceflows + frf_rt, fff_rt, flf_rt = get_structured_faceflows( + flowja, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) + np.testing.assert_array_almost_equal(frf_rt[:, :, :-1], qright[:, :, :-1]) + np.testing.assert_array_almost_equal(fff_rt[:, :-1, :], qfront[:, :-1, :]) def test_get_structured_flowja_to_connections_multilayer(): @@ -826,16 +849,30 @@ def test_get_structured_flowja_to_connections_multilayer(): (qright, qfront, qlower), ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol ) - # first cell of top layer should connect to lower layer + # Cell (0,0,0): lower connection is outflow from n → negative node = 0 # (k=0, i=0, j=0) node_below = nrow * ncol # (k=1, i=0, j=0) conns = ja[ia[node] : ia[node + 1]] flows = flowja[ia[node] : ia[node + 1]] + conn_map = dict(zip(conns, flows)) + assert conn_map[node_below] == -30.0 # lower outflow: -qlower[0,0,0] - # check lower connection - lower_idx = np.where(conns == node_below)[0] - if len(lower_idx) > 0: - assert flows[lower_idx[0]] == 30.0 + # Cell (k=1, i=0, j=0): upper connection is inflow to n from above → positive + node_mid = nrow * ncol # same cell, middle layer + node_above = 0 + conns_mid = ja[ia[node_mid] : ia[node_mid + 1]] + flows_mid = flowja[ia[node_mid] : ia[node_mid + 1]] + conn_map_mid = dict(zip(conns_mid, flows_mid)) + assert conn_map_mid[node_above] == 30.0 # upper inflow: +qlower[0,0,0] + + # Verify roundtrip (interior faces only — boundary cells have no stored neighbor) + from flopy.mf6.utils import get_structured_faceflows + frf_rt, fff_rt, flf_rt = get_structured_faceflows( + flowja, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol + ) + np.testing.assert_array_almost_equal(frf_rt[:, :, :-1], qright[:, :, :-1]) + np.testing.assert_array_almost_equal(fff_rt[:, :-1, :], qfront[:, :-1, :]) + np.testing.assert_array_almost_equal(flf_rt[:-1, :, :], qlower[:-1, :, :]) def test_get_saturation_confined_cells(): @@ -963,7 +1000,8 @@ def test_get_structured_connectivity_big_grid(): # But we only store upper triangle, so expect 4 (diagonal + right + front + lower) node = 1 * nrow * ncol + 10 * ncol + 10 # Middle of layer 1 nconn = ia[node + 1] - ia[node] - assert nconn == 4, f"Interior cell should have 4 connections, got {nconn}" + # Full bidirectional: diagonal + right + left + front + back + lower + upper = 7 + assert nconn == 7, f"Interior cell should have 7 connections, got {nconn}" def test_get_structured_flowja_invalid_inputs(): diff --git a/flopy/mf6/utils/postprocessing.py b/flopy/mf6/utils/postprocessing.py index 47843988fe..2c6cbd8ad5 100644 --- a/flopy/mf6/utils/postprocessing.py +++ b/flopy/mf6/utils/postprocessing.py @@ -33,23 +33,24 @@ def get_structured_connectivity(nlay, nrow, ncol, idomain=None): Notes ----- - Connectivity order for structured grids (upper triangle only): - 1. Diagonal (self connection) - 2. Right (+1 in j, same k, i) - 3. Front (+1 in i, same k, j) - 4. Lower (+1 in k, same i, j) + The IA/JA arrays encode full bidirectional connectivity as MODFLOW 6 + requires: every connection between cells n and m appears twice, once + in cell n's row and once in cell m's row. For each cell, the diagonal + (self-connection) is stored first, followed by all neighbors in + ascending node-number order. The IA/JA arrays use 0-based indexing (Python convention). - When writing to MF6 binary files, add 1 to convert to Fortran 1-based indexing. + When writing to MF6 binary files, add 1 to convert to Fortran 1-based + indexing. Examples -------- >>> import numpy as np - >>> from flopy.mf6.utils import build_structured_connectivity + >>> from flopy.mf6.utils import get_structured_connectivity >>> nlay, nrow, ncol = 2, 3, 3 - >>> ia, ja, nja = build_structured_connectivity(nlay, nrow, ncol) + >>> ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) >>> print(f"Total cells: {nlay * nrow * ncol}, connections: {nja}") - Total cells: 18, connections: 42 + Total cells: 18, connections: 84 """ ncells = nlay * nrow * ncol @@ -64,47 +65,44 @@ def get_structured_connectivity(nlay, nrow, ncol, idomain=None): f"({nlay}, {nrow}, {ncol})" ) - ia = np.zeros(ncells + 1, dtype=np.int32) - ja_list = [] - nja = 0 - + # First pass: collect all neighbor pairs bidirectionally. + # Only iterate right/front/lower to avoid processing each pair twice, + # but register the connection in both cells' adjacency sets. + adjacency = [set() for _ in range(ncells)] for k in range(nlay): for i in range(nrow): for j in range(ncol): - node = k * nrow * ncol + i * ncol + j - - # Skip inactive cells - they still get an entry in IA if idomain[k, i, j] <= 0: - ia[node + 1] = nja continue + n = k * nrow * ncol + i * ncol + j + for dk, di, dj in ((0, 0, 1), (0, 1, 0), (1, 0, 0)): + k2, i2, j2 = k + dk, i + di, j + dj + if k2 < nlay and i2 < nrow and j2 < ncol: + if idomain[k2, i2, j2] > 0: + m = k2 * nrow * ncol + i2 * ncol + j2 + adjacency[n].add(m) + adjacency[m].add(n) + + # Second pass: convert adjacency sets to CSR. + # Each cell's entry begins with the diagonal (self), then neighbors + # in ascending node-number order. + ia = np.zeros(ncells + 1, dtype=np.int32) + ja_list = [] + nja = 0 - # Add diagonal (self connection) - ja_list.append(node) - nja += 1 - - # Add connections to neighbors (upper triangle only) - # Right neighbor (j+1) - if j + 1 < ncol and idomain[k, i, j + 1] > 0: - m = k * nrow * ncol + i * ncol + (j + 1) - ja_list.append(m) - nja += 1 - - # Front neighbor (i+1) - if i + 1 < nrow and idomain[k, i + 1, j] > 0: - m = k * nrow * ncol + (i + 1) * ncol + j - ja_list.append(m) - nja += 1 - - # Lower neighbor (k+1) - if k + 1 < nlay and idomain[k + 1, i, j] > 0: - m = (k + 1) * nrow * ncol + i * ncol + j - ja_list.append(m) - nja += 1 - - ia[node + 1] = nja + for n in range(ncells): + k, i, j = np.unravel_index(n, (nlay, nrow, ncol)) + if idomain[k, i, j] <= 0: + ia[n + 1] = nja + continue + ja_list.append(n) # diagonal + nja += 1 + for m in sorted(adjacency[n]): + ja_list.append(m) + nja += 1 + ia[n + 1] = nja ja = np.array(ja_list, dtype=np.int32) - return ia, ja, nja @@ -350,37 +348,58 @@ def get_structured_flowja( flowja = np.zeros(nja, dtype=np.float64) for n in range(ncells): - # Skip inactive cells k, i, j = np.unravel_index(n, (nlay, nrow, ncol)) if idomain[k, i, j] <= 0: continue - # Get connections for this cell - istart = ia[n] - iend = ia[n + 1] - - for ipos in range(istart, iend): + for ipos in range(ia[n], ia[n + 1]): m = ja[ipos] - # Diagonal - no self flow if m == n: - flowja[ipos] = 0.0 + flowja[ipos] = 0.0 # diagonal (residual; zero here) continue - # Determine connection type by comparing node numbers km, im, jm = np.unravel_index(m, (nlay, nrow, ncol)) - # Right connection (j increases by 1) + # MODFLOW 6 sign convention for FLOW-JA-FACE: + # flowja[n→m] < 0 means flow leaving n (outflow from n to m) + # flowja[n→m] > 0 means flow entering n (inflow to n from m) + # This matches get_structured_faceflows which does frf[n] = -flowja[n→m_right]. + # + # For each face-flow direction the source array stores the flow + # leaving the lower-numbered cell in the positive direction: + # qright[k,i,j] = flow leaving cell n through its right face + # qfront[k,i,j] = flow leaving cell n through its front face + # qlower[k,i,j] = flow leaving cell n through its lower face + # + # Upper-triangle entry (n→m, m is the higher-indexed neighbor): + # flowja = -face_flow_at_n (outflow from n → negative) + # Lower-triangle entry (n→m, m is the lower-indexed neighbor): + # flowja = +face_flow_at_m (inflow to n from m → positive) + + # Right (m is to the right of n): outflow from n if km == k and im == i and jm == j + 1: - flowja[ipos] = qright[k, i, j] + flowja[ipos] = -qright[k, i, j] + + # Left (m is to the left of n): inflow to n from m + elif km == k and im == i and jm == j - 1: + flowja[ipos] = qright[k, i, jm] - # Front connection (i increases by 1) + # Front (m is in front of n): outflow from n elif km == k and im == i + 1 and jm == j: - flowja[ipos] = qfront[k, i, j] + flowja[ipos] = -qfront[k, i, j] - # Lower connection (k increases by 1) + # Back (m is behind n): inflow to n from m + elif km == k and im == i - 1 and jm == j: + flowja[ipos] = qfront[k, im, j] + + # Lower (m is below n): outflow from n elif km == k + 1 and im == i and jm == j: - flowja[ipos] = qlower[k, i, j] + flowja[ipos] = -qlower[k, i, j] + + # Upper (m is above n): inflow to n from m + elif km == k - 1 and im == i and jm == j: + flowja[ipos] = qlower[km, i, j] return flowja diff --git a/flopy/utils/nwt_to_mf6.py b/flopy/utils/nwt_to_mf6.py index 8ac5c0149b..bf8de5982d 100644 --- a/flopy/utils/nwt_to_mf6.py +++ b/flopy/utils/nwt_to_mf6.py @@ -29,29 +29,38 @@ def get_icelltype_from_laytyp(laytyp): Notes ----- - In MODFLOW-NWT, LAYTYP indicates layer properties: + In MODFLOW-NWT (UPW) and MODFLOW-2005 (LPF), LAYTYP indicates layer + properties: - LAYTYP = 0: Confined, transmissivity constant - - LAYTYP > 0: Convertible, transmissivity varies with saturation + - LAYTYP > 0: Convertible with wetting active + - LAYTYP < 0: Convertible without wetting (rewetting disabled) + + Any non-zero LAYTYP means the layer is convertible regardless of sign. + The sign only controls whether the wetting option is active, not whether + the layer can desaturate. (Note: in LPF with THICKSTRT active, negative + LAYTYP layers use starting-head-based thickness rather than the full + cell thickness, but the layer is still treated as convertible in terms + of transmissivity variation with saturation.) In MODFLOW 6, ICELLTYPE similarly indicates cell properties: - ICELLTYPE = 0: Confined - ICELLTYPE > 0: Convertible For conversion, we use a simple mapping: - - LAYTYP = 0 → ICELLTYPE = 0 - - LAYTYP > 0 → ICELLTYPE = 1 + - LAYTYP = 0 → ICELLTYPE = 0 + - LAYTYP != 0 → ICELLTYPE = 1 Examples -------- >>> import numpy as np >>> from flopy.utils.nwt_to_mf6 import get_icelltype_from_laytyp - >>> laytyp = np.array([1, 0, 0]) # Top layer convertible + >>> laytyp = np.array([1, 0, -1]) # Top and bottom layers convertible >>> icelltype = get_icelltype_from_laytyp(laytyp) >>> print(icelltype) - [1 0 0] + [1 0 1] """ laytyp = np.atleast_1d(laytyp).astype(np.int32) - icelltype = np.where(laytyp > 0, 1, 0).astype(np.int32) + icelltype = np.where(laytyp != 0, 1, 0).astype(np.int32) return icelltype @@ -284,20 +293,37 @@ def convert( return {"grb": grb_path, "hds": hds_path, "bud": bud_path} + def _build_header_lookup(self): + """ + Build a dict mapping (kstp, kper) -> header record from the head file. + + Returns + ------- + dict + Keys are (kstp, kper) tuples, values are header records with + fields including 'pertim' and 'totim'. + """ + # get_kstpkper() returns 0-based indices; recordarray stores 1-based + # file values. Subtract 1 so the lookup keys match self.kstpkper. + return { + (int(rec["kstp"]) - 1, int(rec["kper"]) - 1): rec + for rec in self.hds_obj.recordarray + } + def _write_heads(self, filename, precision, verbose): """Write head file with all time steps.""" from .binaryfile import HeadFile - # Read all heads + header_lookup = self._build_header_lookup() head_dict = {} totim_dict = {} pertim_dict = {} for idx, kstpkper in enumerate(self.kstpkper): head = self.hds_obj.get_data(kstpkper=kstpkper) - totim = self.times[idx] - # Assume pertim = totim for now (could be refined) - pertim = totim + rec = header_lookup[kstpkper] + totim = float(rec["totim"]) + pertim = float(rec["pertim"]) head_dict[kstpkper] = head totim_dict[kstpkper] = totim @@ -339,11 +365,22 @@ def _write_budget(self, filename, precision, verbose): # Build list of records records = [] + header_lookup = self._build_header_lookup() + for idx, kstpkper in enumerate(self.kstpkper): kstp, kper = kstpkper - totim = self.times[idx] - pertim = totim # Simplified - delt = 1.0 # Simplified - could calculate from times + rec = header_lookup[kstpkper] + totim = float(rec["totim"]) + pertim = float(rec["pertim"]) + + # delt: for the first time step within a period pertim equals delt; + # for subsequent steps subtract the previous step's pertim. + # kstp is 0-based here (from get_kstpkper()). + if kstp == 0: + delt = pertim + else: + prev_rec = header_lookup.get((kstp - 1, kper)) + delt = pertim - float(prev_rec["pertim"]) if prev_rec else pertim if verbose: print(f" Processing time step {kstpkper}...") From aaf15eceb34e4b8ac76a3a0cfdb01ff2b652e07f Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 15 Apr 2026 19:27:39 -0700 Subject: [PATCH 03/10] should be more general? --- autotest/test_nwt_to_mf6.py | 227 ++++++++++++++++++++++++++++++++++-- flopy/utils/__init__.py | 6 +- flopy/utils/nwt_to_mf6.py | 70 ++++++----- 3 files changed, 262 insertions(+), 41 deletions(-) diff --git a/autotest/test_nwt_to_mf6.py b/autotest/test_nwt_to_mf6.py index 650201d17e..0bde892f93 100644 --- a/autotest/test_nwt_to_mf6.py +++ b/autotest/test_nwt_to_mf6.py @@ -165,7 +165,7 @@ def test_mfgrdfile_write_validation(): def test_nwt_to_mf6_converter_init(function_tmpdir): import numpy as np - from flopy.utils import NwtToMf6Converter + from flopy.utils import ClassicMfToMf6Converter from flopy.utils.binaryfile import CellBudgetFile, HeadFile nlay, nrow, ncol = 2, 3, 4 @@ -229,7 +229,7 @@ def test_nwt_to_mf6_converter_init(function_tmpdir): ncol=ncol, ) - converter = NwtToMf6Converter( + converter = ClassicMfToMf6Converter( str(hds_file), str(cbc_file), nlay, @@ -265,7 +265,7 @@ def test_nwt_to_mf6_converter_convert(function_tmpdir): get_structured_connectivity, get_structured_faceflows, ) - from flopy.utils import NwtToMf6Converter + from flopy.utils import ClassicMfToMf6Converter from flopy.utils.binaryfile import CellBudgetFile, HeadFile nlay, nrow, ncol = 2, 3, 4 @@ -313,7 +313,7 @@ def test_nwt_to_mf6_converter_convert(function_tmpdir): ] CellBudgetFile.write(str(cbc_file), bud_data, nlay=nlay, nrow=nrow, ncol=ncol) - converter = NwtToMf6Converter( + converter = ClassicMfToMf6Converter( str(hds_file), str(cbc_file), nlay, @@ -368,7 +368,7 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): ModflowRch, ModflowUpw, ) - from flopy.utils import NwtToMf6Converter + from flopy.utils import ClassicMfToMf6Converter from flopy.utils.binaryfile import CellBudgetFile, HeadFile modelname = "watertable" @@ -443,7 +443,7 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): assert hds_file.exists(), "Head file not created" assert cbc_file.exists(), "Budget file not created" - converter = NwtToMf6Converter( + converter = ClassicMfToMf6Converter( str(hds_file), str(cbc_file), nlay, @@ -509,7 +509,7 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): ModflowUpw, ModflowWel, ) - from flopy.utils import NwtToMf6Converter + from flopy.utils import ClassicMfToMf6Converter from flopy.utils.binaryfile import CellBudgetFile, HeadFile modelname = "multilayer" @@ -578,7 +578,7 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): assert hds_file.exists() assert cbc_file.exists() - converter = NwtToMf6Converter( + converter = ClassicMfToMf6Converter( str(hds_file), str(cbc_file), nlay, @@ -671,3 +671,214 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): # Verify saturation is correct for convertible layer sat_data = bud_mf6.get_data(text="DATA-SAT", idx=0) assert sat_data is not None + + +# --------------------------------------------------------------------------- +# Fixtures for example-data models +# --------------------------------------------------------------------------- + + +@pytest.fixture +def freyberg_multilayer_path(example_data_path): + return example_data_path / "freyberg_multilayer_transient" + + +@pytest.fixture +def mf2005_freyberg_path(example_data_path): + return example_data_path / "freyberg" + + +# --------------------------------------------------------------------------- +# Integration tests against real pre-computed model output +# --------------------------------------------------------------------------- + + +@pytest.mark.slow +def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_tmpdir): + """ + Convert the pre-computed freyberg_multilayer_transient NWT outputs to MF6 + and verify that the face-flow roundtrip is exact. + + This test requires no executable — it uses the binary output files already + committed to the repo. It checks: + - GRB file has the correct grid dimensions + - Head values are preserved for all time steps + - get_structured_faceflows(flowja_converted) recovers the original NWT + FLOW RIGHT/FRONT/LOWER FACE values (interior faces only) + - DATA-SAT is present for all time steps + """ + from flopy.mf6.utils import MfGrdFile, get_structured_faceflows + from flopy.modflow import Modflow + from flopy.utils import ClassicMfToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + model_ws = freyberg_multilayer_path + + # Load model to get grid parameters (no exe needed) + mf = Modflow.load( + "freyberg.nam", + model_ws=str(model_ws), + check=False, + load_only=["dis", "upw"], + ) + nlay = int(mf.dis.nlay) + nrow = int(mf.dis.nrow) + ncol = int(mf.dis.ncol) + + converter = ClassicMfToMf6Converter( + hds_file=str(model_ws / "freyberg.hds"), + cbc_file=str(model_ws / "freyberg.cbc"), + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=mf.dis.delr.array, + delc=mf.dis.delc.array, + top=mf.dis.top.array, + botm=mf.dis.botm.array, + laytyp=mf.upw.laytyp.array, + hdry=float(mf.upw.hdry), + ) + + output_dir = function_tmpdir / "mf6" + result = converter.convert(str(output_dir)) + + # --- GRB dimensions --- + grb = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + hds_orig = HeadFile(str(model_ws / "freyberg.hds")) + hds_mf6 = HeadFile(str(result["hds"])) + cbc_orig = CellBudgetFile(str(model_ws / "freyberg.cbc")) + bud_mf6 = CellBudgetFile(str(result["bud"])) + + kstpkper_list = hds_orig.get_kstpkper() + + # --- Verify a sample of time steps (first, middle, last) --- + check_indices = [0, len(kstpkper_list) // 2, len(kstpkper_list) - 1] + for check_idx in check_indices: + kstpkper = kstpkper_list[check_idx] + + # Head values are preserved exactly + head_orig = hds_orig.get_data(kstpkper=kstpkper) + head_mf6 = hds_mf6.get_data(kstpkper=kstpkper) + np.testing.assert_array_almost_equal( + head_mf6, head_orig, decimal=5, + err_msg=f"Head mismatch at kstpkper={kstpkper}" + ) + + # Face-flow roundtrip: faceflows → flowja → faceflows should recover + # the original values on interior faces. + frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] + flf_orig = cbc_orig.get_data(text="FLOW LOWER FACE", kstpkper=kstpkper)[0] + + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] + frf_rt, fff_rt, flf_rt = get_structured_faceflows( + flowja, grb_file=str(result["grb"]) + ) + + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5, + err_msg=f"FLOW RIGHT FACE mismatch at kstpkper={kstpkper}" + ) + np.testing.assert_array_almost_equal( + fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5, + err_msg=f"FLOW FRONT FACE mismatch at kstpkper={kstpkper}" + ) + np.testing.assert_array_almost_equal( + flf_rt[:-1, :, :], flf_orig[:-1, :, :], decimal=5, + err_msg=f"FLOW LOWER FACE mismatch at kstpkper={kstpkper}" + ) + + # DATA-SAT is present + sat = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) + assert sat is not None, f"DATA-SAT missing at kstpkper={kstpkper}" + + +@requires_exe("mf2005") +@pytest.mark.slow +def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): + """ + Run the single-layer steady-state Freyberg MODFLOW-2005 model, convert + its outputs to MF6, and verify correctness. + + Tests that the converter works with the LPF package (vs UPW for NWT), + confirming ClassicMfToMf6Converter handles both code paths. + """ + import shutil + + from flopy.mf6.utils import MfGrdFile, get_structured_faceflows + from flopy.modflow import Modflow + from flopy.utils import ClassicMfToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + # Copy model to a writable temp directory and run it + run_ws = function_tmpdir / "mf2005" + shutil.copytree(mf2005_freyberg_path, run_ws) + + mf = Modflow.load( + "freyberg.nam", + model_ws=str(run_ws), + exe_name="mf2005", + check=False, + ) + success, _ = mf.run_model(silent=True) + assert success, "MODFLOW-2005 freyberg run failed" + + hds_path = run_ws / "freyberg.hds" + cbc_path = run_ws / "freyberg.cbc" + assert hds_path.exists(), "Head file not produced" + assert cbc_path.exists(), "Budget file not produced" + + nlay = int(mf.dis.nlay) + nrow = int(mf.dis.nrow) + ncol = int(mf.dis.ncol) + + converter = ClassicMfToMf6Converter( + hds_file=str(hds_path), + cbc_file=str(cbc_path), + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=mf.dis.delr.array, + delc=mf.dis.delc.array, + top=mf.dis.top.array, + botm=mf.dis.botm.array, + laytyp=mf.lpf.laytyp.array, + ) + + output_dir = function_tmpdir / "mf6" + result = converter.convert(str(output_dir)) + + grb = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + hds_orig = HeadFile(str(hds_path)) + hds_mf6 = HeadFile(str(result["hds"])) + cbc_orig = CellBudgetFile(str(cbc_path)) + bud_mf6 = CellBudgetFile(str(result["bud"])) + + kstpkper = hds_orig.get_kstpkper()[0] + + # Head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), + hds_orig.get_data(kstpkper=kstpkper), + decimal=5, + ) + + # Face-flow roundtrip (single layer, so no lower face) + frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] + + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] + frf_rt, fff_rt, _ = get_structured_faceflows( + flowja, grb_file=str(result["grb"]) + ) + + np.testing.assert_array_almost_equal(frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5) + np.testing.assert_array_almost_equal(fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5) diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 5d628fc860..08e94145fb 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -40,7 +40,11 @@ from .mfreadnam import parsenamefile from .modpathfile import EndpointFile, PathlineFile, TimeseriesFile from .mtlistfile import MtListBudget -from .nwt_to_mf6 import NwtToMf6Converter, get_icelltype_from_laytyp +from .nwt_to_mf6 import ( + ClassicMfToMf6Converter, + NwtToMf6Converter, # backward-compatible alias + get_icelltype_from_laytyp, +) from .observationfile import HydmodObs, Mf6Obs, SwrObs from .optionblock import OptionBlock from .postprocessing import get_specific_discharge, get_transmissivities diff --git a/flopy/utils/nwt_to_mf6.py b/flopy/utils/nwt_to_mf6.py index bf8de5982d..4d73a35d33 100644 --- a/flopy/utils/nwt_to_mf6.py +++ b/flopy/utils/nwt_to_mf6.py @@ -64,20 +64,25 @@ def get_icelltype_from_laytyp(laytyp): return icelltype -class NwtToMf6Converter: +class ClassicMfToMf6Converter: """ - Convert MODFLOW-NWT binary outputs to MODFLOW 6 binary format. + Convert classic MODFLOW binary outputs to MODFLOW 6 binary format. - This class reads NWT head and budget files, applies necessary - transformations, and writes MF6-compatible binary files that can - be consumed by PRT or other MF6 post-processors. + Supports any structured-grid MODFLOW variant that produces compact + budget files with FLOW RIGHT FACE / FLOW FRONT FACE / FLOW LOWER FACE + terms — including MODFLOW-NWT (UPW package), MODFLOW-2005 (LPF package), + MODFLOW-2000, and MODFLOW-OWHM. + + This class reads head and cell-budget files, applies the transformations + needed to express flows in MODFLOW 6 format, and writes MF6-compatible + binary files consumable by PRT and other MF6 post-processors. Parameters ---------- hds_file : str - Path to NWT head file (.hds) + Path to head file (.hds) cbc_file : str - Path to NWT cell budget file (.cbc) + Path to cell budget file (.cbc) nlay : int Number of layers nrow : int @@ -93,11 +98,12 @@ class NwtToMf6Converter: botm : array_like Bottom elevation, shape (nlay, nrow, ncol) laytyp : array_like - Layer type from LPF/UPW, shape (nlay,). - 0 = confined, >0 = convertible + Layer type from LPF or UPW package, shape (nlay,). + 0 = confined; any non-zero value = convertible (sign controls + wetting, not confinement). idomain : array_like, optional Domain array, shape (nlay, nrow, ncol). - >0 = active, 0 = inactive, <0 = vertical pass-through + >0 = active, 0 = inactive, <0 = vertical pass-through. If None, all cells are active. hdry : float, optional Head value for dry cells (default -999.0) @@ -109,28 +115,20 @@ class NwtToMf6Converter: Examples -------- >>> import numpy as np - >>> from flopy.utils import NwtToMf6Converter - >>> # Set up grid parameters + >>> from flopy.utils import ClassicMfToMf6Converter >>> nlay, nrow, ncol = 3, 10, 10 >>> delr = np.ones(ncol) * 100.0 >>> delc = np.ones(nrow) * 100.0 >>> top = np.ones((nrow, ncol)) * 10.0 >>> botm = np.zeros((nlay, nrow, ncol)) - >>> botm[0] = 5.0 - >>> botm[1] = 0.0 - >>> botm[2] = -5.0 - >>> laytyp = np.array([1, 0, 0]) # Top layer convertible + >>> botm[0] = 5.0; botm[1] = 0.0; botm[2] = -5.0 + >>> laytyp = np.array([1, 0, 0]) # top layer convertible >>> - >>> # Create converter - >>> converter = NwtToMf6Converter( - ... 'model.hds', - ... 'model.cbc', + >>> converter = ClassicMfToMf6Converter( + ... 'model.hds', 'model.cbc', ... nlay, nrow, ncol, - ... delr, delc, top, botm, - ... laytyp + ... delr, delc, top, botm, laytyp ... ) - >>> - >>> # Convert files >>> converter.convert('mf6_output') """ @@ -325,9 +323,13 @@ def _write_heads(self, filename, precision, verbose): totim = float(rec["totim"]) pertim = float(rec["pertim"]) - head_dict[kstpkper] = head - totim_dict[kstpkper] = totim - pertim_dict[kstpkper] = pertim + # Write 1-based (kstp, kper) so MF6 readers that add 1 internally + # find the correct records. + kstp, kper = kstpkper + kstpkper_out = (int(kstp) + 1, int(kper) + 1) + head_dict[kstpkper_out] = head + totim_dict[kstpkper_out] = totim + pertim_dict[kstpkper_out] = pertim if verbose: print( @@ -477,8 +479,8 @@ def _write_budget(self, filename, precision, verbose): records.append( { "data": flowja, - "kstp": kstp, - "kper": kper, + "kstp": int(kstp) + 1, + "kper": int(kper) + 1, "totim": totim, "pertim": pertim, "delt": delt, @@ -522,8 +524,8 @@ def _write_budget(self, filename, precision, verbose): records.append( { "data": sat_data, - "kstp": kstp, - "kper": kper, + "kstp": int(kstp) + 1, + "kper": int(kper) + 1, "totim": totim, "pertim": pertim, "delt": delt, @@ -549,10 +551,14 @@ def _write_budget(self, filename, precision, verbose): def __repr__(self): return ( - f"NwtToMf6Converter(\n" + f"{type(self).__name__}(\n" f" hds_file={self.hds_file},\n" f" cbc_file={self.cbc_file},\n" f" grid={self.nlay}x{self.nrow}x{self.ncol},\n" f" time_steps={len(self.times)}\n" f")" ) + + +#: Backward-compatible alias. New code should use ClassicMfToMf6Converter. +NwtToMf6Converter = ClassicMfToMf6Converter From dad7d169b55eeee25bfe1f164ffeeee99e5abd31 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 06:21:06 -0700 Subject: [PATCH 04/10] more testing --- autotest/test_nwt_to_mf6.py | 36 +++++++++++++++++++++------------ autotest/test_postprocessing.py | 16 ++++++++------- 2 files changed, 32 insertions(+), 20 deletions(-) diff --git a/autotest/test_nwt_to_mf6.py b/autotest/test_nwt_to_mf6.py index 0bde892f93..50edab07da 100644 --- a/autotest/test_nwt_to_mf6.py +++ b/autotest/test_nwt_to_mf6.py @@ -764,8 +764,10 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t head_orig = hds_orig.get_data(kstpkper=kstpkper) head_mf6 = hds_mf6.get_data(kstpkper=kstpkper) np.testing.assert_array_almost_equal( - head_mf6, head_orig, decimal=5, - err_msg=f"Head mismatch at kstpkper={kstpkper}" + head_mf6, + head_orig, + decimal=5, + err_msg=f"Head mismatch at kstpkper={kstpkper}", ) # Face-flow roundtrip: faceflows → flowja → faceflows should recover @@ -780,16 +782,22 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t ) np.testing.assert_array_almost_equal( - frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5, - err_msg=f"FLOW RIGHT FACE mismatch at kstpkper={kstpkper}" + frf_rt[:, :, :-1], + frf_orig[:, :, :-1], + decimal=5, + err_msg=f"FLOW RIGHT FACE mismatch at kstpkper={kstpkper}", ) np.testing.assert_array_almost_equal( - fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5, - err_msg=f"FLOW FRONT FACE mismatch at kstpkper={kstpkper}" + fff_rt[:, :-1, :], + fff_orig[:, :-1, :], + decimal=5, + err_msg=f"FLOW FRONT FACE mismatch at kstpkper={kstpkper}", ) np.testing.assert_array_almost_equal( - flf_rt[:-1, :, :], flf_orig[:-1, :, :], decimal=5, - err_msg=f"FLOW LOWER FACE mismatch at kstpkper={kstpkper}" + flf_rt[:-1, :, :], + flf_orig[:-1, :, :], + decimal=5, + err_msg=f"FLOW LOWER FACE mismatch at kstpkper={kstpkper}", ) # DATA-SAT is present @@ -876,9 +884,11 @@ def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] - frf_rt, fff_rt, _ = get_structured_faceflows( - flowja, grb_file=str(result["grb"]) - ) + frf_rt, fff_rt, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) - np.testing.assert_array_almost_equal(frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5) - np.testing.assert_array_almost_equal(fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5) + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 + ) + np.testing.assert_array_almost_equal( + fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5 + ) diff --git a/autotest/test_postprocessing.py b/autotest/test_postprocessing.py index 41cb8c5f5a..dfaa98cb4d 100644 --- a/autotest/test_postprocessing.py +++ b/autotest/test_postprocessing.py @@ -813,25 +813,26 @@ def test_get_structured_flowja_to_connections_simple(): conns = ja[ia[node] : ia[node + 1]] flows = flowja[ia[node] : ia[node + 1]] conn_map = dict(zip(conns, flows)) - assert conn_map[node] == 0.0 # diagonal - assert conn_map[1] == -1.0 # right outflow: -qright[0,0,0] - assert conn_map[3] == -2.0 # front outflow: -qfront[0,0,0] + assert conn_map[node] == 0.0 # diagonal + assert conn_map[1] == -1.0 # right outflow: -qright[0,0,0] + assert conn_map[3] == -2.0 # front outflow: -qfront[0,0,0] # Interior cell (0,1,1): left/back/right/front all present node = 1 * ncol + 1 # (k=0, i=1, j=1) conns = ja[ia[node] : ia[node + 1]] flows = flowja[ia[node] : ia[node + 1]] conn_map = dict(zip(conns, flows)) - assert conn_map[node] == 0.0 # diagonal - assert conn_map[node - 1] == qright[0, 1, 0] # left inflow: +qright[0,1,0] + assert conn_map[node] == 0.0 # diagonal + assert conn_map[node - 1] == qright[0, 1, 0] # left inflow: +qright[0,1,0] assert conn_map[node + 1] == -qright[0, 1, 1] # right outflow: -qright[0,1,1] - assert conn_map[node - ncol] == qfront[0, 0, 1] # back inflow: +qfront[0,0,1] - assert conn_map[node + ncol] == -qfront[0, 1, 1] # front outflow: -qfront[0,1,1] + assert conn_map[node - ncol] == qfront[0, 0, 1] # back inflow: +qfront[0,0,1] + assert conn_map[node + ncol] == -qfront[0, 1, 1] # front outflow: -qfront[0,1,1] # Verify roundtrip: flowja → get_structured_faceflows → should recover originals. # Boundary cells with no neighbor in that direction have no stored connection, # so the roundtrip can only recover values for interior faces. from flopy.mf6.utils import get_structured_faceflows + frf_rt, fff_rt, flf_rt = get_structured_faceflows( flowja, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol ) @@ -867,6 +868,7 @@ def test_get_structured_flowja_to_connections_multilayer(): # Verify roundtrip (interior faces only — boundary cells have no stored neighbor) from flopy.mf6.utils import get_structured_faceflows + frf_rt, fff_rt, flf_rt = get_structured_faceflows( flowja, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol ) From 51e4e01b929624a29c634d4dc7fe97abe1d43b1f Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 12:43:14 -0700 Subject: [PATCH 05/10] better naming, explanation of compact files --- ...t_nwt_to_mf6.py => test_classic_to_mf6.py} | 636 ++++++++++++++++-- flopy/utils/__init__.py | 11 +- .../{nwt_to_mf6.py => classic_to_mf6.py} | 346 ++++++++-- 3 files changed, 891 insertions(+), 102 deletions(-) rename autotest/{test_nwt_to_mf6.py => test_classic_to_mf6.py} (55%) rename flopy/utils/{nwt_to_mf6.py => classic_to_mf6.py} (56%) diff --git a/autotest/test_nwt_to_mf6.py b/autotest/test_classic_to_mf6.py similarity index 55% rename from autotest/test_nwt_to_mf6.py rename to autotest/test_classic_to_mf6.py index 50edab07da..40009f41e8 100644 --- a/autotest/test_nwt_to_mf6.py +++ b/autotest/test_classic_to_mf6.py @@ -4,7 +4,7 @@ def test_get_icelltype_from_laytyp(): - from flopy.utils.nwt_to_mf6 import get_icelltype_from_laytyp + from flopy.utils.classic_to_mf6 import get_icelltype_from_laytyp # 1D array laytyp = np.array([1, 0, 0, 1]) @@ -33,6 +33,30 @@ def test_get_icelltype_from_laytyp(): assert np.array_equal(icelltype, [0, 1, 1, 1, 1]) +def test_get_icelltype_from_laycon(): + from flopy.utils.classic_to_mf6 import get_icelltype_from_laycon + + # scalar confined (laycon=0) + assert get_icelltype_from_laycon(0) == 0 + # scalar convertible (laycon=1) + assert get_icelltype_from_laycon(1) == 1 + # laycon=2 is head-dependent-T confined — maps to 0 + assert get_icelltype_from_laycon(2) == 0 + # laycon=3 is fully convertible — maps to 1 + assert get_icelltype_from_laycon(3) == 1 + + # full 4-value array covering all cases + laycon = np.array([0, 1, 2, 3]) + result = get_icelltype_from_laycon(laycon) + assert np.array_equal(result, [0, 1, 0, 1]) + + # shape is preserved + laycon_3d = np.array([[[0, 1], [2, 3]]]) + result_3d = get_icelltype_from_laycon(laycon_3d) + assert result_3d.shape == laycon_3d.shape + assert np.array_equal(result_3d, [[[0, 1], [0, 1]]]) + + def test_mfgrdfile_write_roundtrip(tmp_path): from flopy.mf6.utils import MfGrdFile, get_structured_connectivity @@ -443,19 +467,7 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): assert hds_file.exists(), "Head file not created" assert cbc_file.exists(), "Budget file not created" - converter = ClassicMfToMf6Converter( - str(hds_file), - str(cbc_file), - nlay, - nrow, - ncol, - delr=np.ones(ncol) * delr, - delc=np.ones(nrow) * delc, - top=np.ones((nrow, ncol)) * top, - botm=np.ones((nlay, nrow, ncol)) * botm, - laytyp=np.array([1]), # Convertible - model_ws=function_tmpdir, - ) + converter = ClassicMfToMf6Converter.from_model(mf, hds_file, cbc_file) output_dir = function_tmpdir / "mf6_output" result = converter.convert(str(output_dir), verbose=True) @@ -489,10 +501,24 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): assert "DATA-SAT" in texts_str, "DATA-SAT not in budget" flowja = bud_mf6.get_data(text="FLOW-JA-FACE", idx=0) - assert flowja is not None, "Could not read FLOW-JA-FACE" - - sat_data = bud_mf6.get_data(text="DATA-SAT", idx=0) - assert sat_data is not None, "Could not read DATA-SAT" + assert len(flowja) > 0, "Could not read FLOW-JA-FACE" + + kstpkper0 = bud_mf6.get_kstpkper()[0] + sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper0) + assert len(sat_data) > 0, "Could not read DATA-SAT" + + # Validate saturation values numerically. For the single convertible layer + # (laytyp=1, icelltype=1) with top=25 and botm=0: + # sat = clamp((head - botm) / (top - botm), 0, 1) = clamp(head / 25, 0, 1) + # DATA-SAT is stored as imeth=6; the saturation values are in the "q" field. + sat_vals = sat_data[0]["q"] + expected_sat = np.clip(head_nwt.flatten() / (top - botm), 0.0, 1.0) + np.testing.assert_array_almost_equal( + sat_vals, + expected_sat, + decimal=5, + err_msg="DATA-SAT values do not match expected (head - botm) / (top - botm)", + ) @requires_exe("mfnwt") @@ -578,19 +604,7 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): assert hds_file.exists() assert cbc_file.exists() - converter = ClassicMfToMf6Converter( - str(hds_file), - str(cbc_file), - nlay, - nrow, - ncol, - delr=np.ones(ncol) * delr, - delc=np.ones(nrow) * delc, - top=top, - botm=botm, - laytyp=laytyp, - model_ws=function_tmpdir, - ) + converter = ClassicMfToMf6Converter.from_model(mf, hds_file, cbc_file) output_dir = function_tmpdir / "mf6_output" result = converter.convert(str(output_dir), verbose=True) @@ -669,8 +683,9 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): assert flowja is not None # Verify saturation is correct for convertible layer - sat_data = bud_mf6.get_data(text="DATA-SAT", idx=0) - assert sat_data is not None + kstpkper0 = bud_mf6.get_kstpkper()[0] + sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper0) + assert len(sat_data) > 0, "Could not read DATA-SAT" # --------------------------------------------------------------------------- @@ -702,10 +717,10 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t This test requires no executable — it uses the binary output files already committed to the repo. It checks: - GRB file has the correct grid dimensions - - Head values are preserved for all time steps + - Head values are preserved for a sample of time steps (first, middle, last) - get_structured_faceflows(flowja_converted) recovers the original NWT FLOW RIGHT/FRONT/LOWER FACE values (interior faces only) - - DATA-SAT is present for all time steps + - DATA-SAT is present for sampled time steps """ from flopy.mf6.utils import MfGrdFile, get_structured_faceflows from flopy.modflow import Modflow @@ -802,7 +817,7 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t # DATA-SAT is present sat = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) - assert sat is not None, f"DATA-SAT missing at kstpkper={kstpkper}" + assert len(sat) > 0, f"DATA-SAT missing at kstpkper={kstpkper}" @requires_exe("mf2005") @@ -844,6 +859,241 @@ def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): nrow = int(mf.dis.nrow) ncol = int(mf.dis.ncol) + converter = ClassicMfToMf6Converter.from_model(mf, hds_path, cbc_path) + + output_dir = function_tmpdir / "mf6" + result = converter.convert(str(output_dir)) + + grb = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + hds_orig = HeadFile(str(hds_path)) + hds_mf6 = HeadFile(str(result["hds"])) + cbc_orig = CellBudgetFile(str(cbc_path)) + bud_mf6 = CellBudgetFile(str(result["bud"])) + + kstpkper = hds_orig.get_kstpkper()[0] + + # Head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), + hds_orig.get_data(kstpkper=kstpkper), + decimal=5, + ) + + # Face-flow roundtrip (single layer, so no lower face) + frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] + + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] + frf_rt, fff_rt, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) + + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 + ) + np.testing.assert_array_almost_equal( + fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5 + ) + + +@requires_exe("mf2000") +@pytest.mark.slow +def test_classic_to_mf6_mf2000_watertable(function_tmpdir): + """ + Build and run a 1-layer unconfined watertable MODFLOW-2000 model, convert + its outputs to MF6 format, and verify head, face-flow, and saturation. + + MODFLOW-2000 produces compact budget files with the same binary format as + MODFLOW-NWT and MODFLOW-2005. This test confirms ClassicMfToMf6Converter + handles MF-2000 output identically to later variants. + """ + from flopy.mf6.utils import MfGrdFile, get_structured_faceflows + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowGhb, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowRch, + ) + from flopy.utils import ClassicMfToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + modelname = "mf2000_watertable" + nlay, nrow, ncol = 1, 1, 100 + delr, delc = 50.0, 1.0 + h1, h2 = 20.0, 11.0 + top, botm = 25.0, 0.0 + hk = 50.0 + + # Linear initial heads to avoid convergence failures with PCG and unconfined cells + strt = np.zeros((nlay, nrow, ncol)) + strt[0, 0, :] = np.linspace(h1, h2, ncol) + + # GHB boundary conditions at both ends + h_adj1 = h1 - (h1 - h2) / ncol + h_adj2 = h2 + (h1 - h2) / ncol + b1 = 0.5 * (h1 + h_adj1) + b2 = 0.5 * (h2 + h_adj2) + c1 = hk * b1 * delc / (0.5 * delr) + c2 = hk * b2 * delc / (0.5 * delr) + ghb_dtype = ModflowGhb.get_default_dtype() + spd = np.zeros(2, dtype=ghb_dtype).view(np.recarray) + spd[0] = (0, 0, 0, h1, c1) + spd[1] = (0, 0, ncol - 1, h2, c2) + + mf = Modflow( + modelname=modelname, + exe_name="mf2000", + model_ws=str(function_tmpdir), + version="mf2k", + ) + ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + perlen=1, + nstp=1, + steady=True, + ) + ModflowBas(mf, ibound=1, strt=strt) + ModflowLpf(mf, hk=hk, laytyp=1, ipakcb=53) + ModflowGhb(mf, stress_period_data=spd) + ModflowRch(mf, rech=0.001, nrchop=1) + oc = ModflowOc( + mf, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, + ) + oc.reset_budgetunit(budgetunit=53, fname=f"{modelname}.cbc") + ModflowPcg(mf) + + mf.write_input() + success, _ = mf.run_model(silent=True) + assert success, "MODFLOW-2000 watertable run failed" + + hds_path = function_tmpdir / f"{modelname}.hds" + cbc_path = function_tmpdir / f"{modelname}.cbc" + assert hds_path.exists() + assert cbc_path.exists() + + converter = ClassicMfToMf6Converter.from_model(mf, hds_path, cbc_path) + result = converter.convert(str(function_tmpdir / "mf6_output")) + + assert result["grb"].exists() + assert result["hds"].exists() + assert result["bud"].exists() + + hds_orig = HeadFile(str(hds_path)) + hds_mf6 = HeadFile(str(result["hds"])) + cbc_orig = CellBudgetFile(str(cbc_path)) + bud_mf6 = CellBudgetFile(str(result["bud"])) + + kstpkper = hds_orig.get_kstpkper()[0] + + # Head roundtrip + head_orig = hds_orig.get_data(kstpkper=kstpkper) + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 + ) + + # Face-flow roundtrip (1 row, 1 layer → only right face flows exist) + frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] + frf_rt, _, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 + ) + + # Saturation: sat = clamp(head / (top - botm), 0, 1) + sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) + assert len(sat_data) > 0, "DATA-SAT missing" + expected_sat = np.clip(head_orig.flatten() / (top - botm), 0.0, 1.0) + np.testing.assert_array_almost_equal(sat_data[0]["q"], expected_sat, decimal=5) + + +@requires_exe("mf2005") +@pytest.mark.slow +@pytest.mark.parametrize( + "model_key", + [ + "mf6/test/test001a_Tharmonic", # 1 layer, 1 stress period, steady + "mf6/test/test001h_rch_array3", # 1 layer, 4 stress periods + ], +) +def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): + """ + Fetch a MODFLOW-2005/LPF companion model from the devtools registry, + run it, convert outputs to MF6, and verify head and face-flow roundtrip. + + The mf6/test registry includes mf2005/ subdirectories alongside each MF6 + test case. Running these models exercises the LPF code path across a + variety of boundary conditions and confirms multi-stress-period conversion. + """ + import shutil + + from modflow_devtools.models import copy_to + + from flopy.mf6.utils import get_structured_faceflows + from flopy.modflow import Modflow + from flopy.utils import ClassicMfToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + # Download the parent MF6 model; the mf2005/ companion is a subdirectory. + model_ws = function_tmpdir / "model" + copy_to(str(model_ws), model_key) + mf2005_ws = model_ws / "mf2005" + assert mf2005_ws.exists(), f"mf2005/ subdir not found for {model_key}" + + # Locate the namefile (there should be exactly one .nam in mf2005/) + nam_files = list(mf2005_ws.glob("*.nam")) + assert len(nam_files) == 1, f"Expected 1 .nam file in mf2005/, got {nam_files}" + nam_name = nam_files[0].name + + # Run mf2005 + mf = Modflow.load( + nam_name, + model_ws=str(mf2005_ws), + exe_name="mf2005", + check=False, + ) + + # The companion models may not have "save budget" in their OC. Force it so + # that the CBC file is populated for the converter. + nper = int(mf.dis.nper) + nstp = mf.dis.nstp.array + sp_data = { + (kstp, kper): ["save head", "save budget"] + for kper in range(nper) + for kstp in range(nstp[kper]) + } + from flopy.modflow import ModflowOc + + ModflowOc(mf, stress_period_data=sp_data, compact=True) + mf.write_input() + + success, _ = mf.run_model(silent=True) + assert success, f"mf2005 run failed for {model_key}" + + # Locate binary outputs written by the OC package + hds_path = next(mf2005_ws.glob("*.hds"), None) + cbc_path = next(mf2005_ws.glob("*.cbc"), None) + assert hds_path is not None, "No .hds file produced" + assert cbc_path is not None, "No .cbc file produced" + + nlay = int(mf.dis.nlay) + nrow = int(mf.dis.nrow) + ncol = int(mf.dis.ncol) + converter = ClassicMfToMf6Converter( hds_file=str(hds_path), cbc_file=str(cbc_path), @@ -856,39 +1106,335 @@ def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): botm=mf.dis.botm.array, laytyp=mf.lpf.laytyp.array, ) + result = converter.convert(str(function_tmpdir / "mf6_output")) - output_dir = function_tmpdir / "mf6" - result = converter.convert(str(output_dir)) + assert result["grb"].exists() + assert result["hds"].exists() + assert result["bud"].exists() + + hds_orig = HeadFile(str(hds_path)) + hds_mf6 = HeadFile(str(result["hds"])) + cbc_orig = CellBudgetFile(str(cbc_path)) + bud_mf6 = CellBudgetFile(str(result["bud"])) + + for kstpkper in hds_orig.get_kstpkper(): + # Head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), + hds_orig.get_data(kstpkper=kstpkper), + decimal=5, + err_msg=f"Head mismatch at kstpkper={kstpkper}", + ) + + # Face-flow roundtrip (check only terms that the model actually wrote) + cbc_texts = [ + t.decode().strip() if isinstance(t, bytes) else t.strip() + for t in cbc_orig.textlist + ] + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] + frf_rt, fff_rt, flf_rt = get_structured_faceflows( + flowja, grb_file=str(result["grb"]) + ) + if "FLOW RIGHT FACE" in cbc_texts: + frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], + frf_orig[:, :, :-1], + decimal=5, + err_msg=f"FLOW RIGHT FACE mismatch at kstpkper={kstpkper}", + ) + if "FLOW FRONT FACE" in cbc_texts: + fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] + np.testing.assert_array_almost_equal( + fff_rt[:, :-1, :], + fff_orig[:, :-1, :], + decimal=5, + err_msg=f"FLOW FRONT FACE mismatch at kstpkper={kstpkper}", + ) + if "FLOW LOWER FACE" in cbc_texts: + flf_orig = cbc_orig.get_data(text="FLOW LOWER FACE", kstpkper=kstpkper)[0] + np.testing.assert_array_almost_equal( + flf_rt[:-1, :, :], + flf_orig[:-1, :, :], + decimal=5, + err_msg=f"FLOW LOWER FACE mismatch at kstpkper={kstpkper}", + ) + + # DATA-SAT present + sat = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) + assert len(sat) > 0, f"DATA-SAT missing at kstpkper={kstpkper}" + + +@requires_exe("mf2005") +@pytest.mark.slow +def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): + """ + Build and run a 3-layer MODFLOW-2005/LPF model (1 convertible + 2 confined), + convert its outputs to MF6, and verify head, face-flow, and saturation. + + Complements test_classic_to_mf6_freyberg_mf2005 (single layer) and the + NWT multilayer tests by exercising the LPF code path with nlay > 1. + """ + from flopy.mf6.utils import MfGrdFile, get_structured_faceflows + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowRch, + ModflowWel, + ) + from flopy.utils import ClassicMfToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + modelname = "mf2005_multilayer" + nlay, nrow, ncol = 3, 10, 10 + delr = delc = 100.0 + top = np.ones((nrow, ncol)) * 100.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 80.0 + botm[1] = 60.0 + botm[2] = 40.0 + laytyp = np.array([1, 0, 0]) # top convertible, lower two confined + strt = np.ones((nlay, nrow, ncol)) * 95.0 + hk = np.ones((nlay, nrow, ncol)) * 10.0 + + mf = Modflow( + modelname=modelname, + exe_name="mf2005", + model_ws=str(function_tmpdir), + version="mf2005", + ) + ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + perlen=1, + nstp=1, + steady=True, + ) + ModflowBas(mf, ibound=1, strt=strt) + ModflowLpf(mf, hk=hk, laytyp=laytyp, ipakcb=53) + ModflowWel(mf, stress_period_data={0: [(1, 5, 5, -1000.0)]}) + ModflowRch(mf, rech=0.001) + oc = ModflowOc( + mf, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, + ) + oc.reset_budgetunit(budgetunit=53, fname=f"{modelname}.cbc") + ModflowPcg(mf) + + mf.write_input() + success, _ = mf.run_model(silent=True) + assert success, "MODFLOW-2005 multilayer run failed" + + hds_path = function_tmpdir / f"{modelname}.hds" + cbc_path = function_tmpdir / f"{modelname}.cbc" + assert hds_path.exists() + assert cbc_path.exists() + + converter = ClassicMfToMf6Converter.from_model(mf, hds_path, cbc_path) + result = converter.convert(str(function_tmpdir / "mf6_output")) + + # GRB dimensions grb = MfGrdFile(str(result["grb"])) assert grb.nlay == nlay assert grb.nrow == nrow assert grb.ncol == ncol + # ICELLTYPE pattern (Fortran interleaved: layer varies fastest) + icelltype = grb._datadict["ICELLTYPE"] + assert np.all(icelltype[0::nlay] == 1), "Layer 1 should be convertible" + assert np.all(icelltype[1::nlay] == 0), "Layer 2 should be confined" + assert np.all(icelltype[2::nlay] == 0), "Layer 3 should be confined" + hds_orig = HeadFile(str(hds_path)) hds_mf6 = HeadFile(str(result["hds"])) cbc_orig = CellBudgetFile(str(cbc_path)) bud_mf6 = CellBudgetFile(str(result["bud"])) kstpkper = hds_orig.get_kstpkper()[0] + head_orig = hds_orig.get_data(kstpkper=kstpkper) # Head roundtrip np.testing.assert_array_almost_equal( - hds_mf6.get_data(kstpkper=kstpkper), - hds_orig.get_data(kstpkper=kstpkper), - decimal=5, + hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 ) - # Face-flow roundtrip (single layer, so no lower face) + # Face-flow roundtrip (3 layers → include lower face) frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] - + flf_orig = cbc_orig.get_data(text="FLOW LOWER FACE", kstpkper=kstpkper)[0] flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] - frf_rt, fff_rt, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) - + frf_rt, fff_rt, flf_rt = get_structured_faceflows( + flowja, grb_file=str(result["grb"]) + ) np.testing.assert_array_almost_equal( frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 ) np.testing.assert_array_almost_equal( fff_rt[:, :-1, :], fff_orig[:, :-1, :], decimal=5 ) + np.testing.assert_array_almost_equal( + flf_rt[:-1, :, :], flf_orig[:-1, :, :], decimal=5 + ) + + # Saturation: confined layers = 1.0; convertible layer = (head - botm) / thickness + sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) + assert len(sat_data) > 0, "DATA-SAT missing" + sat_flat = sat_data[0]["q"] # Fortran order: layer varies fastest + # Confined layers (1 and 2, 0-indexed) should all be 1.0 + np.testing.assert_array_equal( + sat_flat[1::nlay], 1.0, err_msg="Layer 2 (confined) sat != 1" + ) + np.testing.assert_array_equal( + sat_flat[2::nlay], 1.0, err_msg="Layer 3 (confined) sat != 1" + ) + # Convertible layer sat should match (head - botm[0]) / (top - botm[0]) + head_layer0 = head_orig[0].flatten(order="F") + top_layer0 = top.flatten(order="F") + bot_layer0 = botm[0].flatten(order="F") + expected_sat_layer0 = np.clip( + (head_layer0 - bot_layer0) / (top_layer0 - bot_layer0), 0.0, 1.0 + ) + np.testing.assert_array_almost_equal( + sat_flat[0::nlay], expected_sat_layer0, decimal=5 + ) + + +@requires_exe("mf2000") +@pytest.mark.slow +def test_classic_to_mf6_mf2000_bcf(function_tmpdir): + """ + Build and run a MODFLOW-2000 model with the BCF package (laycon=3, + fully convertible), convert via from_model(), and verify that + get_icelltype_from_laycon() correctly maps laycon → ICELLTYPE and that + head, face-flow, and saturation all round-trip correctly. + + This test exercises the BCF code path in from_model(), which differs + from the LPF/UPW path: laycon values (0/1/2/3) are first converted to + ICELLTYPE (0/1) via get_icelltype_from_laycon(), and hdry is read from + the BCF package (default -1e30, not -999.0). + """ + from flopy.mf6.utils import get_structured_faceflows + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowBcf, + ModflowDis, + ModflowGhb, + ModflowOc, + ModflowPcg, + ModflowRch, + ) + from flopy.utils import ClassicMfToMf6Converter + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + modelname = "mf2000_bcf" + nlay, nrow, ncol = 1, 1, 100 + delr, delc = 50.0, 1.0 + h1, h2 = 20.0, 11.0 + top, botm = 25.0, 0.0 + hk = 50.0 + + # Linear initial heads to avoid convergence failures with PCG and unconfined cells + strt = np.zeros((nlay, nrow, ncol)) + strt[0, 0, :] = np.linspace(h1, h2, ncol) + + h_adj1 = h1 - (h1 - h2) / ncol + h_adj2 = h2 + (h1 - h2) / ncol + b1, b2 = 0.5 * (h1 + h_adj1), 0.5 * (h2 + h_adj2) + c1 = hk * b1 * delc / (0.5 * delr) + c2 = hk * b2 * delc / (0.5 * delr) + ghb_dtype = ModflowGhb.get_default_dtype() + spd = np.zeros(2, dtype=ghb_dtype).view(np.recarray) + spd[0] = (0, 0, 0, h1, c1) + spd[1] = (0, 0, ncol - 1, h2, c2) + + mf = Modflow( + modelname=modelname, + exe_name="mf2000", + model_ws=str(function_tmpdir), + version="mf2k", + ) + ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + perlen=1, + nstp=1, + steady=True, + ) + ModflowBas(mf, ibound=1, strt=strt) + # laycon=3: fully convertible (T varies with saturated thickness) + ModflowBcf(mf, laycon=3, hy=hk, ipakcb=53) + ModflowGhb(mf, stress_period_data=spd) + ModflowRch(mf, rech=0.001, nrchop=1) + oc = ModflowOc( + mf, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, + ) + oc.reset_budgetunit(budgetunit=53, fname=f"{modelname}.cbc") + ModflowPcg(mf) + + mf.write_input() + success, _ = mf.run_model(silent=True) + assert success, "MODFLOW-2000 BCF run failed" + + hds_path = function_tmpdir / f"{modelname}.hds" + cbc_path = function_tmpdir / f"{modelname}.cbc" + assert hds_path.exists() + assert cbc_path.exists() + + # from_model() must detect BCF, call get_icelltype_from_laycon(laycon=3)→1, + # and read hdry from bcf6 (default -1e30 for BCF). + converter = ClassicMfToMf6Converter.from_model(mf, hds_path, cbc_path) + assert converter.icelltype[0] == 1, "laycon=3 should give icelltype=1" + assert converter.hdry == pytest.approx(-1e30), "hdry should come from BCF package" + + result = converter.convert(str(function_tmpdir / "mf6_output")) + assert result["grb"].exists() + assert result["hds"].exists() + assert result["bud"].exists() + + hds_orig = HeadFile(str(hds_path)) + hds_mf6 = HeadFile(str(result["hds"])) + cbc_orig = CellBudgetFile(str(cbc_path)) + bud_mf6 = CellBudgetFile(str(result["bud"])) + + kstpkper = hds_orig.get_kstpkper()[0] + head_orig = hds_orig.get_data(kstpkper=kstpkper) + + # Head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 + ) + + # Face-flow roundtrip (1 row, 1 layer → only right face flows exist) + frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] + frf_rt, _, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 + ) + + # Saturation (laycon=3 → convertible → sat = clamp(head / 25, 0, 1)) + sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) + assert len(sat_data) > 0, "DATA-SAT missing" + active = head_orig.flatten() > -1e29 # exclude hdry/hnoflo + expected_sat = np.clip(head_orig.flatten()[active] / (top - botm), 0.0, 1.0) + np.testing.assert_array_almost_equal(sat_data[0]["q"], expected_sat, decimal=5) diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 08e94145fb..3026f43bb4 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -29,6 +29,12 @@ from .formattedfile import FormattedHeadFile get_modflow = get_modflow_module.run_main +from .classic_to_mf6 import ( + ClassicMfToMf6Converter, + NwtToMf6Converter, # backward-compatible alias + get_icelltype_from_laycon, + get_icelltype_from_laytyp, +) from .gridintersect import GridIntersect from .mflistfile import ( Mf6ListBudget, @@ -40,11 +46,6 @@ from .mfreadnam import parsenamefile from .modpathfile import EndpointFile, PathlineFile, TimeseriesFile from .mtlistfile import MtListBudget -from .nwt_to_mf6 import ( - ClassicMfToMf6Converter, - NwtToMf6Converter, # backward-compatible alias - get_icelltype_from_laytyp, -) from .observationfile import HydmodObs, Mf6Obs, SwrObs from .optionblock import OptionBlock from .postprocessing import get_specific_discharge, get_transmissivities diff --git a/flopy/utils/nwt_to_mf6.py b/flopy/utils/classic_to_mf6.py similarity index 56% rename from flopy/utils/nwt_to_mf6.py rename to flopy/utils/classic_to_mf6.py index 4d73a35d33..e2b3e770d1 100644 --- a/flopy/utils/nwt_to_mf6.py +++ b/flopy/utils/classic_to_mf6.py @@ -1,5 +1,69 @@ """ -Utilities for converting MODFLOW-NWT binary outputs to MODFLOW 6 format. +Utilities for converting classic MODFLOW binary outputs to MODFLOW 6 format. + +Supported variants +------------------ +Any structured-grid MODFLOW variant that writes compact budget files with +FLOW RIGHT FACE / FLOW FRONT FACE / FLOW LOWER FACE terms is supported, +including MODFLOW-NWT (UPW), MODFLOW-2005 (LPF), and MODFLOW-2000 (LPF or +BCF). Use :class:`ClassicMfToMf6Converter` or its +:meth:`~ClassicMfToMf6Converter.from_model` +constructor to perform the conversion. + +Binary format differences: classic MODFLOW vs MODFLOW 6 +-------------------------------------------------------- +Understanding what this converter does requires knowing how the two file +formats differ. + +**Head files (.hds)** + The on-disk layout is essentially identical: one record per layer per + time step, each prefixed by a header containing ``(kstp, kper, pertim, + totim, text, ncol, nrow, ilay)``. The only practical difference is that + MF6 readers add 1 to the caller-supplied ``(kstp, kper)`` before + searching, so the values written in the file must be 1-based. Classic + MODFLOW writes 1-based values natively, but flopy's ``get_kstpkper()`` + returns 0-based pairs; this converter re-adds 1 before writing. + +**Budget files (.cbc / .bud)** + This is where the formats diverge significantly. + + Classic MODFLOW budget files exist in two sub-formats controlled by the + OC package keyword ``COMPACT BUDGET [FILES]``: + + * *Non-compact* (old-style, ``COMPACT BUDGET`` absent): ``nlay > 0`` in + every record header; no extended header fields; all records are + implicitly full 3-D arrays (``imeth=0``); no timing information is + stored. **This converter does not support the non-compact format.** + + * *Compact* (``COMPACT BUDGET`` or ``COMPACT BUDGET FILES`` in OC): + ``nlay < 0`` is the binary signal; each record has an extended header + with ``(imeth, delt, pertim, totim)``; face-flow terms use ``imeth=1`` + (full 3-D array) and boundary-flux terms use ``imeth=2`` or higher. + **This is the format the converter requires.** + + Within the compact format, classic and MF6 semantics diverge: + + * *Classic compact* stores each flux component as a separate named + record—``FLOW RIGHT FACE``, ``FLOW FRONT FACE``, ``FLOW LOWER FACE``— + each a 3-D array of shape ``(nlay, nrow, ncol)`` (``imeth=1``). + Inter-cell flows are directional and half-edge (each face stored once, + positive away from the lower-index cell). + + * *MF6* stores all inter-cell flows in a single ``FLOW-JA-FACE`` record: + a 1-D array of length NJA (total number of cell connections) indexed by + the IA/JA sparse connectivity arrays. Each connection appears twice + (once from each side) with opposite signs. Saturation and specific + discharge are stored as separate ``DATA-SAT`` and ``DATA-SPDIS`` + records using the ``imeth=6`` sparse-list format. MF6 also uses + ``nlay < 0``, so it shares the same compact record framing—just + entirely different content. + +**Grid file (.grb) — new in MF6** + Classic MODFLOW stores grid geometry in the ASCII DIS file. MF6 + requires a binary grid record (``.dis.grb``) that embeds the IA/JA + connectivity, cell geometry (TOP, BOT, DELR, DELC), ICELLTYPE, and + IDOMAIN. MF6 post-processors (PRT, etc.) use the GRB to interpret + ``FLOW-JA-FACE``; there is no classic-MODFLOW equivalent. """ import os @@ -10,12 +74,12 @@ def get_icelltype_from_laytyp(laytyp): """ - Convert NWT LAYTYP values to MF6 ICELLTYPE values. + Convert classic MODFLOW LAYTYP values to MF6 ICELLTYPE values. Parameters ---------- laytyp : array_like - Layer type array from NWT LPF or UPW package. + Layer type array from a UPW (NWT) or LPF (2005/2000) package. - 0: Confined - >0: Convertible (unconfined/water table) Shape can be (nlay,) or (nlay, nrow, ncol) @@ -29,7 +93,7 @@ def get_icelltype_from_laytyp(laytyp): Notes ----- - In MODFLOW-NWT (UPW) and MODFLOW-2005 (LPF), LAYTYP indicates layer + In MODFLOW-NWT (UPW) and MODFLOW-2005/2000 (LPF), LAYTYP indicates layer properties: - LAYTYP = 0: Confined, transmissivity constant - LAYTYP > 0: Convertible with wetting active @@ -53,7 +117,7 @@ def get_icelltype_from_laytyp(laytyp): Examples -------- >>> import numpy as np - >>> from flopy.utils.nwt_to_mf6 import get_icelltype_from_laytyp + >>> from flopy.utils import get_icelltype_from_laytyp >>> laytyp = np.array([1, 0, -1]) # Top and bottom layers convertible >>> icelltype = get_icelltype_from_laytyp(laytyp) >>> print(icelltype) @@ -64,70 +128,147 @@ def get_icelltype_from_laytyp(laytyp): return icelltype +def get_icelltype_from_laycon(laycon): + """ + Convert BCF LAYCON values to MF6 ICELLTYPE values. + + Parameters + ---------- + laycon : array_like + Layer connectivity array from the BCF package, shape ``(nlay,)`` or + ``(nlay, nrow, ncol)``. + + * 0 — confined; transmissivity is constant. + * 1 — convertible; transmissivity based on initial saturated thickness. + * 2 — confined; transmissivity varies with head (quasi-3D / Tran + specified). + * 3 — convertible; transmissivity varies with saturated thickness. + + Returns + ------- + icelltype : ndarray + Cell type for MF6 (same shape as input): + + * 0 — confined (LAYCON 0 or 2). + * 1 — convertible (LAYCON 1 or 3). + + Notes + ----- + LAYCON 1 and LAYCON 3 both represent convertible layers in which the + water table can fall below the cell top; they differ only in how + transmissivity is computed from the saturated thickness. Both map to + ``ICELLTYPE = 1``. + + LAYCON 2 layers have a head-dependent transmissivity but are never + allowed to desaturate fully, so they map to ``ICELLTYPE = 0`` + (confined). + + Examples + -------- + >>> import numpy as np + >>> from flopy.utils import get_icelltype_from_laycon + >>> laycon = np.array([3, 0, 2, 1]) + >>> get_icelltype_from_laycon(laycon) + array([1, 0, 0, 1], dtype=int32) + """ + laycon = np.atleast_1d(laycon).astype(np.int32) + icelltype = np.where((laycon == 1) | (laycon == 3), 1, 0).astype(np.int32) + return icelltype + + class ClassicMfToMf6Converter: """ Convert classic MODFLOW binary outputs to MODFLOW 6 binary format. - Supports any structured-grid MODFLOW variant that produces compact - budget files with FLOW RIGHT FACE / FLOW FRONT FACE / FLOW LOWER FACE - terms — including MODFLOW-NWT (UPW package), MODFLOW-2005 (LPF package), - MODFLOW-2000, and MODFLOW-OWHM. + Reads head and cell-budget files produced by any structured-grid MODFLOW + variant with compact budget output (FLOW RIGHT FACE / FLOW FRONT FACE / + FLOW LOWER FACE terms), and writes the MF6-compatible head file, budget + file (``FLOW-JA-FACE``, ``DATA-SAT``), and binary grid record (GRB) + required by PRT and other MF6 post-processors. - This class reads head and cell-budget files, applies the transformations - needed to express flows in MODFLOW 6 format, and writes MF6-compatible - binary files consumable by PRT and other MF6 post-processors. + Confirmed compatible variants: MODFLOW-NWT (UPW), MODFLOW-2005 (LPF), + MODFLOW-2000 (LPF or BCF). + + The easiest way to construct a converter from a loaded model is + :meth:`from_model`, which autodetects the flow package and extracts all + required parameters automatically. Construct directly only when a model + object is not available (e.g. you only have the binary files on disk). Parameters ---------- - hds_file : str - Path to head file (.hds) - cbc_file : str - Path to cell budget file (.cbc) + hds_file : str or PathLike + Path to the head file (.hds). + cbc_file : str or PathLike + Path to the compact cell-budget file (.cbc). nlay : int - Number of layers + Number of layers. nrow : int - Number of rows + Number of rows. ncol : int - Number of columns - delr : array_like - Column spacing, shape (ncol,) - delc : array_like - Row spacing, shape (nrow,) - top : array_like - Top elevation, shape (nrow, ncol) - botm : array_like - Bottom elevation, shape (nlay, nrow, ncol) - laytyp : array_like - Layer type from LPF or UPW package, shape (nlay,). - 0 = confined; any non-zero value = convertible (sign controls - wetting, not confinement). - idomain : array_like, optional - Domain array, shape (nlay, nrow, ncol). - >0 = active, 0 = inactive, <0 = vertical pass-through. - If None, all cells are active. + Number of columns. + delr : array_like, shape (ncol,) + Column widths (along rows). + delc : array_like, shape (nrow,) + Row widths (along columns). + top : array_like, shape (nrow, ncol) + Top elevation of the model. + botm : array_like, shape (nlay, nrow, ncol) + Bottom elevation of each layer. + laytyp : array_like, shape (nlay,) + Layer type from the LPF or UPW package. ``0`` → confined; any + non-zero value → convertible (sign controls wetting, not + confinement). For BCF models, pass the output of + :func:`get_icelltype_from_laycon` here rather than the raw + ``laycon`` values. + idomain : array_like, shape (nlay, nrow, ncol), optional + Active-cell mask: ``> 0`` active, ``0`` inactive, + ``< 0`` vertical pass-through. ``None`` treats all cells as active. hdry : float, optional - Head value for dry cells (default -999.0) + Sentinel head value written by MODFLOW for dry cells. Should match + the value in the LPF/UPW ``HDRY`` field (default ``-999.0``). For + BCF models the default is ``-1e30``; pass the model value explicitly + or use :meth:`from_model`. hnoflo : float, optional - Head value for inactive cells (default -9999.0) + Sentinel head value for inactive (IBOUND ≤ 0) cells. Should match + the BAS6 ``HNOFLO`` value (default ``-9999.0``; BAS6 default is + ``-999.99``). Use :meth:`from_model` to avoid mismatches. model_ws : str or PathLike, optional - Model workspace for input files (default current directory) + Directory prepended to ``hds_file`` and ``cbc_file`` when those + paths are relative. Default is the current directory. + + See Also + -------- + from_model : Preferred constructor when a loaded ``Modflow`` object is + available; handles flow-package detection and sentinel-value + extraction automatically. + get_icelltype_from_laytyp : Maps LPF/UPW ``LAYTYP`` → ``ICELLTYPE``. + get_icelltype_from_laycon : Maps BCF ``LAYCON`` → ``ICELLTYPE``. Examples -------- + Construct from a loaded model (recommended): + + >>> from flopy.modflow import Modflow + >>> from flopy.utils import ClassicMfToMf6Converter + >>> mf = Modflow.load('model.nam', load_only=['dis', 'upw', 'bas6']) + >>> converter = ClassicMfToMf6Converter.from_model( + ... mf, 'model.hds', 'model.cbc' + ... ) + >>> converter.convert('mf6_output') + + Construct directly from arrays (when model object is unavailable): + >>> import numpy as np >>> from flopy.utils import ClassicMfToMf6Converter >>> nlay, nrow, ncol = 3, 10, 10 - >>> delr = np.ones(ncol) * 100.0 - >>> delc = np.ones(nrow) * 100.0 - >>> top = np.ones((nrow, ncol)) * 10.0 - >>> botm = np.zeros((nlay, nrow, ncol)) - >>> botm[0] = 5.0; botm[1] = 0.0; botm[2] = -5.0 >>> laytyp = np.array([1, 0, 0]) # top layer convertible - >>> >>> converter = ClassicMfToMf6Converter( ... 'model.hds', 'model.cbc', ... nlay, nrow, ncol, - ... delr, delc, top, botm, laytyp + ... np.ones(ncol) * 100.0, np.ones(nrow) * 100.0, + ... np.ones((nrow, ncol)) * 10.0, + ... np.array([5.0, 0.0, -5.0])[:, None, None] * np.ones((nlay, nrow, ncol)), + ... laytyp, ... ) >>> converter.convert('mf6_output') """ @@ -202,6 +343,108 @@ def __init__( self.times = self.hds_obj.get_times() self.kstpkper = self.hds_obj.get_kstpkper() + @classmethod + def from_model(cls, model, hds_file, cbc_file, **kwargs): + """ + Construct a converter from a loaded ``Modflow`` model object. + + Autodetects the flow package (UPW, LPF, or BCF), extracts grid + geometry and layer-type arrays, and reads ``hdry`` / ``hnoflo`` + sentinel values so they do not need to be provided manually. + + Parameters + ---------- + model : flopy.modflow.Modflow + Loaded MODFLOW model. At minimum the DIS package and one of + UPW, LPF, or BCF must be loaded. BAS6 is optional but + recommended so that ``hnoflo`` is read from the model rather + than using the default. + hds_file : str or PathLike + Path to the head file produced by the model run. + cbc_file : str or PathLike + Path to the compact cell-budget file produced by the model run. + **kwargs + Additional keyword arguments forwarded to + :class:`ClassicMfToMf6Converter` (e.g. ``idomain``, + ``model_ws``). Parameters already extracted from ``model`` + (``nlay``, ``nrow``, ``ncol``, ``delr``, ``delc``, ``top``, + ``botm``, ``laytyp``, ``hdry``, ``hnoflo``) must not be + supplied here. + + Returns + ------- + ClassicMfToMf6Converter + + Raises + ------ + ValueError + If the model has no DIS package, or has no supported flow + package (UPW, LPF, or BCF). + + Examples + -------- + NWT / UPW model: + + >>> from flopy.modflow import Modflow + >>> from flopy.utils import ClassicMfToMf6Converter + >>> mf = Modflow.load('model.nam', load_only=['dis', 'upw', 'bas6']) + >>> converter = ClassicMfToMf6Converter.from_model( + ... mf, 'model.hds', 'model.cbc' + ... ) + + MF2000 / BCF model: + + >>> mf = Modflow.load('model.nam', load_only=['dis', 'bcf6', 'bas6']) + >>> converter = ClassicMfToMf6Converter.from_model( + ... mf, 'model.hds', 'model.cbc' + ... ) + """ + if model.dis is None: + raise ValueError( + "model.dis is None — load the DIS package before calling from_model()." + ) + + # Detect flow package and derive laytyp / hdry. + # For BCF, laycon values are first converted to ICELLTYPE (0/1) + # and passed as laytyp; __init__ then maps them through + # get_icelltype_from_laytyp, which is a no-op for 0/1 values. + if getattr(model, "upw", None) is not None: + laytyp = model.upw.laytyp.array + hdry = float(model.upw.hdry) + elif getattr(model, "lpf", None) is not None: + laytyp = model.lpf.laytyp.array + hdry = float(model.lpf.hdry) + elif getattr(model, "bcf6", None) is not None: + laytyp = get_icelltype_from_laycon(model.bcf6.laycon.array) + hdry = float(model.bcf6.hdry) + else: + raise ValueError( + "No supported flow package found on model. " + "Load one of: UPW, LPF, or BCF before calling from_model()." + ) + + # Read hnoflo from BAS6 if available. + if getattr(model, "bas6", None) is not None: + hnoflo = float(model.bas6.hnoflo) + else: + hnoflo = -9999.0 + + return cls( + hds_file=str(hds_file), + cbc_file=str(cbc_file), + nlay=int(model.dis.nlay), + nrow=int(model.dis.nrow), + ncol=int(model.dis.ncol), + delr=model.dis.delr.array, + delc=model.dis.delc.array, + top=model.dis.top.array, + botm=model.dis.botm.array, + laytyp=laytyp, + hdry=hdry, + hnoflo=hnoflo, + **kwargs, + ) + def convert( self, output_dir, @@ -212,12 +455,12 @@ def convert( verbose=False, ): """ - Convert NWT binary outputs to MF6 format. + Convert classic MODFLOW binary outputs to MF6 format. Parameters ---------- output_dir : str or PathLike - Directory for output files (will be created if doesn't exist) + Directory for output files (will be created if it does not exist) grb_name : str, optional Grid file name (default 'gwf.grb') hds_name : str, optional @@ -249,7 +492,7 @@ def convert( bud_path = output_dir / bud_name if verbose: - print("\nConverting NWT outputs to MF6 format") + print("\nConverting classic MODFLOW outputs to MF6 format") print(f" Input head file: {self.hds_file}") print(f" Input budget file: {self.cbc_file}") print(f" Output directory: {output_dir}") @@ -490,10 +733,9 @@ def _write_budget(self, filename, precision, verbose): ) # 2. DATA-SPDIS (specific discharge) - SKIPPED for now - # TODO: Requires refactoring get_specific_discharge() to work - # without model object or implementing a minimal wrapper. PRT can - # calculate specific discharge from FLOW-JA-FACE if needed. - # See PHASE3_REMAINING_ISSUES.md for details. + # TODO: get_specific_discharge() requires a model object and cannot + # easily be driven from binary files alone. PRT can reconstruct + # specific discharge from FLOW-JA-FACE if needed. # 3. Calculate saturation sat = get_saturation( From 43cf18606e3626f1be58f684f549805d54cba947 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 14:52:04 -0700 Subject: [PATCH 06/10] improve comments/docstrings, add tests, just copy head file thru --- autotest/test_classic_to_mf6.py | 93 +++++++++++++++++++++++++++++++++ flopy/utils/classic_to_mf6.py | 91 +++++++------------------------- 2 files changed, 112 insertions(+), 72 deletions(-) diff --git a/autotest/test_classic_to_mf6.py b/autotest/test_classic_to_mf6.py index 40009f41e8..b6feabd658 100644 --- a/autotest/test_classic_to_mf6.py +++ b/autotest/test_classic_to_mf6.py @@ -3,6 +3,99 @@ from modflow_devtools.markers import requires_exe +def test_headfile_classic_readable_without_conversion(tmp_path): + """ + A classic MODFLOW head file (1-based kstpkper in the binary header) can be + read directly by HeadFile without any conversion step. The values returned + by get_data() on the original file must be identical to those returned after + round-tripping through _write_heads. + """ + from flopy.utils.binaryfile import HeadFile + + nlay, nrow, ncol = 2, 3, 4 + rng = np.random.default_rng(0) + head1 = rng.standard_normal((nlay, nrow, ncol)) + head2 = rng.standard_normal((nlay, nrow, ncol)) + + # Write with 1-based kstpkper at known timing, as classic MODFLOW does. + original = tmp_path / "classic.hds" + HeadFile.write( + str(original), + {(1, 1): head1, (2, 1): head2}, + nlay=nlay, + nrow=nrow, + ncol=ncol, + precision="double", + totim={(1, 1): 1.0, (2, 1): 2.0}, + pertim={(1, 1): 1.0, (2, 1): 1.0}, + ) + + hds = HeadFile(str(original)) + # get_kstpkper() returns 0-based — (0,0) and (1,0) for the two records above + kstpkper_list = hds.get_kstpkper() + assert kstpkper_list == [(0, 0), (1, 0)] + + # Data is retrievable by 0-based index — no conversion needed + np.testing.assert_array_equal(hds.get_data(kstpkper=(0, 0)), head1) + np.testing.assert_array_equal(hds.get_data(kstpkper=(1, 0)), head2) + + +def test_headfile_convert_is_bytewise_noop(tmp_path): + """ + convert() copies the head file verbatim — the output .hds is byte-for-byte + identical to the original classic MODFLOW .hds. + + This confirms that the converter does not need to rewrite head files. + """ + from flopy.utils.binaryfile import CellBudgetFile, HeadFile + + from flopy.utils.classic_to_mf6 import ClassicMfToMf6Converter + + nlay, nrow, ncol = 2, 3, 4 + rng = np.random.default_rng(1) + head = rng.standard_normal((nlay, nrow, ncol)) + + original = tmp_path / "classic.hds" + HeadFile.write( + str(original), + {(1, 1): head}, + nlay=nlay, + nrow=nrow, + ncol=ncol, + precision="double", + totim={(1, 1): 1.0}, + pertim={(1, 1): 1.0}, + ) + + delr = np.ones(ncol) * 100.0 + delc = np.ones(nrow) * 100.0 + top = np.ones((nrow, ncol)) * 10.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 5.0 + laytyp = np.array([1, 0]) + + cbc_file = tmp_path / "dummy.cbc" + frf = np.zeros((nlay, nrow, ncol)) + CellBudgetFile.write( + str(cbc_file), + [{"data": frf.flatten(order="F"), "kstp": 1, "kper": 1, + "text": "FLOW RIGHT FACE", "imeth": 1}], + nlay=nlay, nrow=nrow, ncol=ncol, + ) + + converter = ClassicMfToMf6Converter( + str(original), str(cbc_file), + nlay, nrow, ncol, delr, delc, top, botm, laytyp, + ) + + output_dir = tmp_path / "mf6_out" + result = converter.convert(str(output_dir)) + + assert result["hds"].read_bytes() == original.read_bytes(), ( + "convert() changed the head file — expected a verbatim copy" + ) + + def test_get_icelltype_from_laytyp(): from flopy.utils.classic_to_mf6 import get_icelltype_from_laytyp diff --git a/flopy/utils/classic_to_mf6.py b/flopy/utils/classic_to_mf6.py index e2b3e770d1..bae8f29f30 100644 --- a/flopy/utils/classic_to_mf6.py +++ b/flopy/utils/classic_to_mf6.py @@ -1,72 +1,60 @@ """ -Utilities for converting classic MODFLOW binary outputs to MODFLOW 6 format. +Utilities for converting legacy MODFLOW binary outputs to MODFLOW 6 format. Supported variants ------------------ Any structured-grid MODFLOW variant that writes compact budget files with -FLOW RIGHT FACE / FLOW FRONT FACE / FLOW LOWER FACE terms is supported, +FLOW RIGHT FACE, FLOW FRONT FACE, and FLOW LOWER FACE terms is supported, including MODFLOW-NWT (UPW), MODFLOW-2005 (LPF), and MODFLOW-2000 (LPF or BCF). Use :class:`ClassicMfToMf6Converter` or its :meth:`~ClassicMfToMf6Converter.from_model` constructor to perform the conversion. -Binary format differences: classic MODFLOW vs MODFLOW 6 +Binary format differences: compact MODFLOW vs MODFLOW 6 -------------------------------------------------------- -Understanding what this converter does requires knowing how the two file -formats differ. - **Head files (.hds)** The on-disk layout is essentially identical: one record per layer per time step, each prefixed by a header containing ``(kstp, kper, pertim, - totim, text, ncol, nrow, ilay)``. The only practical difference is that - MF6 readers add 1 to the caller-supplied ``(kstp, kper)`` before - searching, so the values written in the file must be 1-based. Classic - MODFLOW writes 1-based values natively, but flopy's ``get_kstpkper()`` - returns 0-based pairs; this converter re-adds 1 before writing. + totim, text, ncol, nrow, ilay)``. **Budget files (.cbc / .bud)** - This is where the formats diverge significantly. - Classic MODFLOW budget files exist in two sub-formats controlled by the OC package keyword ``COMPACT BUDGET [FILES]``: * *Non-compact* (old-style, ``COMPACT BUDGET`` absent): ``nlay > 0`` in every record header; no extended header fields; all records are - implicitly full 3-D arrays (``imeth=0``); no timing information is - stored. **This converter does not support the non-compact format.** + implicitly full 3-D arrays (``imeth=0``); stress period and time step + are stored but no time data. **Not supported by this converter.** * *Compact* (``COMPACT BUDGET`` or ``COMPACT BUDGET FILES`` in OC): ``nlay < 0`` is the binary signal; each record has an extended header with ``(imeth, delt, pertim, totim)``; face-flow terms use ``imeth=1`` - (full 3-D array) and boundary-flux terms use ``imeth=2`` or higher. + (full 3-D array) and boundary-flow terms use ``imeth=2`` or higher. **This is the format the converter requires.** Within the compact format, classic and MF6 semantics diverge: - * *Classic compact* stores each flux component as a separate named + * *Classic compact* stores each flow component as a separate named record—``FLOW RIGHT FACE``, ``FLOW FRONT FACE``, ``FLOW LOWER FACE``— each a 3-D array of shape ``(nlay, nrow, ncol)`` (``imeth=1``). - Inter-cell flows are directional and half-edge (each face stored once, - positive away from the lower-index cell). + Inter-cell flows are directional with each face stored once, in the + positive direction (away from the lower-index cell). * *MF6* stores all inter-cell flows in a single ``FLOW-JA-FACE`` record: a 1-D array of length NJA (total number of cell connections) indexed by the IA/JA sparse connectivity arrays. Each connection appears twice (once from each side) with opposite signs. Saturation and specific discharge are stored as separate ``DATA-SAT`` and ``DATA-SPDIS`` - records using the ``imeth=6`` sparse-list format. MF6 also uses - ``nlay < 0``, so it shares the same compact record framing—just - entirely different content. + records using the ``imeth=6`` sparse-list format. MF6 budget file + headers also use ``nlay < 0``, like classic compact files. **Grid file (.grb) — new in MF6** Classic MODFLOW stores grid geometry in the ASCII DIS file. MF6 - requires a binary grid record (``.dis.grb``) that embeds the IA/JA - connectivity, cell geometry (TOP, BOT, DELR, DELC), ICELLTYPE, and - IDOMAIN. MF6 post-processors (PRT, etc.) use the GRB to interpret - ``FLOW-JA-FACE``; there is no classic-MODFLOW equivalent. + flow models create a binary grid file (``.dis.grb``) that encodes + the IA/JA connectivity, cell geometry, convertibility, and IDOMAIN. """ -import os +import shutil from pathlib import Path import numpy as np @@ -519,10 +507,11 @@ def convert( verbose=verbose, ) - # Write HDS file + # Copy head file verbatim — the on-disk format is identical between + # classic MODFLOW and MF6, so no conversion is needed. if verbose: - print(f"\nWriting head file: {hds_path}") - self._write_heads(hds_path, precision, verbose) + print(f"\nCopying head file to: {hds_path}") + shutil.copy2(self.hds_file, hds_path) # Write BUD file if verbose: @@ -551,48 +540,6 @@ def _build_header_lookup(self): for rec in self.hds_obj.recordarray } - def _write_heads(self, filename, precision, verbose): - """Write head file with all time steps.""" - from .binaryfile import HeadFile - - header_lookup = self._build_header_lookup() - head_dict = {} - totim_dict = {} - pertim_dict = {} - - for idx, kstpkper in enumerate(self.kstpkper): - head = self.hds_obj.get_data(kstpkper=kstpkper) - rec = header_lookup[kstpkper] - totim = float(rec["totim"]) - pertim = float(rec["pertim"]) - - # Write 1-based (kstp, kper) so MF6 readers that add 1 internally - # find the correct records. - kstp, kper = kstpkper - kstpkper_out = (int(kstp) + 1, int(kper) + 1) - head_dict[kstpkper_out] = head - totim_dict[kstpkper_out] = totim - pertim_dict[kstpkper_out] = pertim - - if verbose: - print( - f" Read head for time step {kstpkper}: " - f"min={head.min():.2f}, max={head.max():.2f}" - ) - - # Write using HeadFile.write() - HeadFile.write( - str(filename), - head_dict, - nlay=self.nlay, - nrow=self.nrow, - ncol=self.ncol, - precision=precision, - totim=totim_dict, - pertim=pertim_dict, - verbose=verbose, - ) - def _write_budget(self, filename, precision, verbose): """Write budget file with FLOW-JA-FACE, DATA-SPDIS, and DATA-SAT.""" from ..mf6.utils import get_structured_flowja From e52c423dc1e03b9037776e95ff3c401259c9634a Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 14:57:23 -0700 Subject: [PATCH 07/10] ruff --- autotest/test_classic_to_mf6.py | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/autotest/test_classic_to_mf6.py b/autotest/test_classic_to_mf6.py index b6feabd658..19a68c064f 100644 --- a/autotest/test_classic_to_mf6.py +++ b/autotest/test_classic_to_mf6.py @@ -48,7 +48,6 @@ def test_headfile_convert_is_bytewise_noop(tmp_path): This confirms that the converter does not need to rewrite head files. """ from flopy.utils.binaryfile import CellBudgetFile, HeadFile - from flopy.utils.classic_to_mf6 import ClassicMfToMf6Converter nlay, nrow, ncol = 2, 3, 4 @@ -78,14 +77,31 @@ def test_headfile_convert_is_bytewise_noop(tmp_path): frf = np.zeros((nlay, nrow, ncol)) CellBudgetFile.write( str(cbc_file), - [{"data": frf.flatten(order="F"), "kstp": 1, "kper": 1, - "text": "FLOW RIGHT FACE", "imeth": 1}], - nlay=nlay, nrow=nrow, ncol=ncol, + [ + { + "data": frf.flatten(order="F"), + "kstp": 1, + "kper": 1, + "text": "FLOW RIGHT FACE", + "imeth": 1, + } + ], + nlay=nlay, + nrow=nrow, + ncol=ncol, ) converter = ClassicMfToMf6Converter( - str(original), str(cbc_file), - nlay, nrow, ncol, delr, delc, top, botm, laytyp, + str(original), + str(cbc_file), + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + laytyp, ) output_dir = tmp_path / "mf6_out" From 0f318f52d67c78353fa3e9d8e73b38abe72b84fc Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 15:19:07 -0700 Subject: [PATCH 08/10] clean up tests --- autotest/test_binaryfile.py | 28 +- autotest/test_binarygrid_util.py | 151 +++++++-- autotest/test_classic_to_mf6.py | 521 ++++--------------------------- 3 files changed, 227 insertions(+), 473 deletions(-) diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index e3308ef9c7..77e75b5956 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -1302,10 +1302,36 @@ def test_headfile_write_scalar_disv_disu(function_tmpdir): hds.close() -def test_empty_list_error(function_tmpdir): +def test_write_empty_list_error(function_tmpdir): """Test that empty lists raise appropriate errors.""" with pytest.raises(ValueError, match="Empty data list"): HeadFile.write(function_tmpdir / "test.hds", []) with pytest.raises(ValueError, match="Empty data list"): CellBudgetFile.write(function_tmpdir / "test.cbc", [], text="STORAGE") + + +def test_headfile_write_roundtrip(function_tmpdir): + nlay, nrow, ncol = 2, 3, 4 + rng = np.random.default_rng(0) + head1 = rng.standard_normal((nlay, nrow, ncol)) + head2 = rng.standard_normal((nlay, nrow, ncol)) + + original = function_tmpdir / "classic.hds" + HeadFile.write( + str(original), + {(1, 1): head1, (2, 1): head2}, + nlay=nlay, + nrow=nrow, + ncol=ncol, + precision="double", + totim={(1, 1): 1.0, (2, 1): 2.0}, + pertim={(1, 1): 1.0, (2, 1): 1.0}, + ) + + hds = HeadFile(str(original)) + kstpkper_list = hds.get_kstpkper() + assert kstpkper_list == [(0, 0), (1, 0)] + + np.testing.assert_array_equal(hds.get_data(kstpkper=(0, 0)), head1) + np.testing.assert_array_equal(hds.get_data(kstpkper=(1, 0)), head2) diff --git a/autotest/test_binarygrid_util.py b/autotest/test_binarygrid_util.py index 5d56e03203..314d5bce5d 100644 --- a/autotest/test_binarygrid_util.py +++ b/autotest/test_binarygrid_util.py @@ -5,7 +5,7 @@ from matplotlib import pyplot as plt from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid -from flopy.mf6.utils import MfGrdFile +from flopy.mf6.utils import MfGrdFile, get_structured_connectivity pytestmark = pytest.mark.mf6 @@ -161,7 +161,7 @@ def test_mfgrddisu_modelgrid(mfgrd_test_path): ) -def test_write_grb_instance_method(tmp_path, mfgrd_test_path): +def test_mfgrdfile_export_instance_method(tmp_path, mfgrd_test_path): original_file = mfgrd_test_path / "nwtp3.dis.grb" grb_orig = MfGrdFile(original_file, verbose=False) @@ -191,7 +191,9 @@ def test_write_grb_instance_method(tmp_path, mfgrd_test_path): np.testing.assert_array_equal(grb_new.idomain, grb_orig.idomain) -def test_write_grb_instance_method_precision_conversion(tmp_path, mfgrd_test_path): +def test_mfgrdfile_export_instance_method_precision_conversion( + tmp_path, mfgrd_test_path +): original_file = mfgrd_test_path / "nwtp3.dis.grb" grb = MfGrdFile(original_file, verbose=False) @@ -209,10 +211,7 @@ def test_write_grb_instance_method_precision_conversion(tmp_path, mfgrd_test_pat assert single_file.stat().st_size < double_file.stat().st_size -def test_write_grb_disv_roundtrip(tmp_path, mfgrd_test_path): - """Test MfGrdFile.export() for DISV grid with roundtrip validation.""" - from flopy.mf6.utils.binarygrid_util import MfGrdFile - +def test_mfgrdfile_export_disv_roundtrip(tmp_path, mfgrd_test_path): # Read original DISV grb file original_file = mfgrd_test_path / "flow.disv.grb" grb_orig = MfGrdFile(original_file, verbose=False) @@ -265,10 +264,7 @@ def test_write_grb_disv_roundtrip(tmp_path, mfgrd_test_path): ) -def test_write_grb_disv_precision_conversion(tmp_path, mfgrd_test_path): - """Test MfGrdFile.export() for DISV grid with precision conversion.""" - from flopy.mf6.utils.binarygrid_util import MfGrdFile - +def test_mfgrdfile_export_disv_precision_conversion(tmp_path, mfgrd_test_path): # Read original DISV grb file original_file = mfgrd_test_path / "flow.disv.grb" grb = MfGrdFile(original_file, verbose=False) @@ -315,10 +311,7 @@ def test_write_grb_disv_precision_conversion(tmp_path, mfgrd_test_path): ) -def test_write_grb_disu_roundtrip(tmp_path, mfgrd_test_path): - """Test MfGrdFile.export() for DISU grid with roundtrip validation.""" - from flopy.mf6.utils.binarygrid_util import MfGrdFile - +def test_mfgrdfile_export_disu_roundtrip(tmp_path, mfgrd_test_path): # Read original DISU grb file original_file = mfgrd_test_path / "flow.disu.grb" grb_orig = MfGrdFile(original_file, verbose=False) @@ -360,10 +353,7 @@ def test_write_grb_disu_roundtrip(tmp_path, mfgrd_test_path): np.testing.assert_array_equal(grb_new.idomain, grb_orig.idomain) -def test_write_grb_disu_precision_conversion(tmp_path, mfgrd_test_path): - """Test MfGrdFile.export() for DISU grid with precision conversion.""" - from flopy.mf6.utils.binarygrid_util import MfGrdFile - +def test_mfgrdfile_export_disu_precision_conversion(tmp_path, mfgrd_test_path): # Read original DISU grb file original_file = mfgrd_test_path / "flow.disu.grb" grb = MfGrdFile(original_file, verbose=False) @@ -396,3 +386,126 @@ def test_write_grb_disu_precision_conversion(tmp_path, mfgrd_test_path): # Double precision should match exactly (same precision as original) np.testing.assert_allclose(grb_double.top, grb.top, rtol=1e-12) np.testing.assert_allclose(grb_double.bot, grb.bot, rtol=1e-12) + + +def test_mfgrdfile_write_dis_roundtrip(tmp_path): + nlay, nrow, ncol = 2, 3, 4 + delr = np.ones(ncol) * 100.0 + delc = np.ones(nrow) * 50.0 + top = np.ones((nrow, ncol)) * 10.0 + botm = np.zeros((nlay, nrow, ncol)) + botm[0] = 5.0 + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) + icelltype[0] = 1 # top layer convertible + + grb_file = tmp_path / "test.dis.grb" + MfGrdFile.write_dis( + str(grb_file), + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + ia, + ja, + icelltype=icelltype, + ) + + grb = MfGrdFile(str(grb_file)) + + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + assert grb.nodes == nlay * nrow * ncol + assert grb.nja == nja + + np.testing.assert_array_almost_equal(grb.delr, delr) + np.testing.assert_array_almost_equal(grb.delc, delc) + + # Build expected TOP for all cells (MF6 format) + # Layer 1: model top, Layer 2+: bottom of layer above + # In Fortran order (layer-interleaved): [L0[0,0], L1[0,0], L2[0,0], L0[0,1], ...] + top_expected = np.zeros(nlay * nrow * ncol) + top_flat = top.flatten(order="F") + botm_flat = botm.flatten(order="F") + + for i in range(nrow * ncol): + # Layer 0: use model top + top_expected[i * nlay] = top_flat[i] + # Layers 1+: use bottom of layer above + for k in range(1, nlay): + top_expected[i * nlay + k] = botm_flat[i * nlay + (k - 1)] + + np.testing.assert_array_almost_equal(grb.top, top_expected) + np.testing.assert_array_almost_equal(grb.bot, botm.flatten(order="F")) + + np.testing.assert_array_equal(grb.ia, ia) + np.testing.assert_array_equal(grb.ja, ja) + + icelltype_read = grb._datadict["ICELLTYPE"] + np.testing.assert_array_equal(icelltype_read, icelltype.flatten(order="F")) + + +def test_mfgrdfile_write_dis_with_idomain(tmp_path): + nlay, nrow, ncol = 1, 3, 3 + delr = np.ones(ncol) + delc = np.ones(nrow) + top = np.ones((nrow, ncol)) + botm = np.zeros((nlay, nrow, ncol)) + + # center cell inactive + idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) + idomain[0, 1, 1] = 0 + + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol, idomain) + + grb_file = tmp_path / "test_idomain.dis.grb" + MfGrdFile.write_dis( + str(grb_file), nlay, nrow, ncol, delr, delc, top, botm, ia, ja, idomain=idomain + ) + + grb = MfGrdFile(str(grb_file)) + idomain_read = grb._datadict["IDOMAIN"] + np.testing.assert_array_equal(idomain_read, idomain.flatten(order="F")) + + +def test_mfgrdfile_write_dis_validation(): + nlay, nrow, ncol = 1, 2, 2 + delr = np.ones(ncol) + delc = np.ones(nrow) + top = np.ones((nrow, ncol)) + botm = np.zeros((nlay, nrow, ncol)) + ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + + with pytest.raises(ValueError, match="delr length"): + MfGrdFile.write_dis( + "test.grb", + nlay, + nrow, + ncol, + np.ones(3), # Wrong length + delc, + top, + botm, + ia, + ja, + ) + + with pytest.raises(ValueError, match="ia length"): + MfGrdFile.write_dis( + "test.grb", + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + np.ones(10), # Wrong length (should be ncells + 1 = 5) + ja, + ) diff --git a/autotest/test_classic_to_mf6.py b/autotest/test_classic_to_mf6.py index 19a68c064f..519f19879e 100644 --- a/autotest/test_classic_to_mf6.py +++ b/autotest/test_classic_to_mf6.py @@ -1,119 +1,24 @@ +import shutil + import numpy as np import pytest from modflow_devtools.markers import requires_exe - -def test_headfile_classic_readable_without_conversion(tmp_path): - """ - A classic MODFLOW head file (1-based kstpkper in the binary header) can be - read directly by HeadFile without any conversion step. The values returned - by get_data() on the original file must be identical to those returned after - round-tripping through _write_heads. - """ - from flopy.utils.binaryfile import HeadFile - - nlay, nrow, ncol = 2, 3, 4 - rng = np.random.default_rng(0) - head1 = rng.standard_normal((nlay, nrow, ncol)) - head2 = rng.standard_normal((nlay, nrow, ncol)) - - # Write with 1-based kstpkper at known timing, as classic MODFLOW does. - original = tmp_path / "classic.hds" - HeadFile.write( - str(original), - {(1, 1): head1, (2, 1): head2}, - nlay=nlay, - nrow=nrow, - ncol=ncol, - precision="double", - totim={(1, 1): 1.0, (2, 1): 2.0}, - pertim={(1, 1): 1.0, (2, 1): 1.0}, - ) - - hds = HeadFile(str(original)) - # get_kstpkper() returns 0-based — (0,0) and (1,0) for the two records above - kstpkper_list = hds.get_kstpkper() - assert kstpkper_list == [(0, 0), (1, 0)] - - # Data is retrievable by 0-based index — no conversion needed - np.testing.assert_array_equal(hds.get_data(kstpkper=(0, 0)), head1) - np.testing.assert_array_equal(hds.get_data(kstpkper=(1, 0)), head2) - - -def test_headfile_convert_is_bytewise_noop(tmp_path): - """ - convert() copies the head file verbatim — the output .hds is byte-for-byte - identical to the original classic MODFLOW .hds. - - This confirms that the converter does not need to rewrite head files. - """ - from flopy.utils.binaryfile import CellBudgetFile, HeadFile - from flopy.utils.classic_to_mf6 import ClassicMfToMf6Converter - - nlay, nrow, ncol = 2, 3, 4 - rng = np.random.default_rng(1) - head = rng.standard_normal((nlay, nrow, ncol)) - - original = tmp_path / "classic.hds" - HeadFile.write( - str(original), - {(1, 1): head}, - nlay=nlay, - nrow=nrow, - ncol=ncol, - precision="double", - totim={(1, 1): 1.0}, - pertim={(1, 1): 1.0}, - ) - - delr = np.ones(ncol) * 100.0 - delc = np.ones(nrow) * 100.0 - top = np.ones((nrow, ncol)) * 10.0 - botm = np.zeros((nlay, nrow, ncol)) - botm[0] = 5.0 - laytyp = np.array([1, 0]) - - cbc_file = tmp_path / "dummy.cbc" - frf = np.zeros((nlay, nrow, ncol)) - CellBudgetFile.write( - str(cbc_file), - [ - { - "data": frf.flatten(order="F"), - "kstp": 1, - "kper": 1, - "text": "FLOW RIGHT FACE", - "imeth": 1, - } - ], - nlay=nlay, - nrow=nrow, - ncol=ncol, - ) - - converter = ClassicMfToMf6Converter( - str(original), - str(cbc_file), - nlay, - nrow, - ncol, - delr, - delc, - top, - botm, - laytyp, - ) - - output_dir = tmp_path / "mf6_out" - result = converter.convert(str(output_dir)) - - assert result["hds"].read_bytes() == original.read_bytes(), ( - "convert() changed the head file — expected a verbatim copy" - ) +from flopy.mf6.utils import ( + MfGrdFile, + get_structured_connectivity, + get_structured_faceflows, +) +from flopy.modflow import Modflow +from flopy.utils.binaryfile import CellBudgetFile, HeadFile +from flopy.utils.classic_to_mf6 import ( + ClassicMfToMf6Converter, + get_icelltype_from_laycon, + get_icelltype_from_laytyp, +) def test_get_icelltype_from_laytyp(): - from flopy.utils.classic_to_mf6 import get_icelltype_from_laytyp # 1D array laytyp = np.array([1, 0, 0, 1]) @@ -138,269 +43,33 @@ def test_get_icelltype_from_laytyp(): # various LAYTYP values — any non-zero means convertible laytyp = np.array([0, 1, 2, 3, -1]) icelltype = get_icelltype_from_laytyp(laytyp) - # 0 stays 0, all others (positive or negative) map to 1 + # 0 stays 0, all others (positive or negative) 1 assert np.array_equal(icelltype, [0, 1, 1, 1, 1]) def test_get_icelltype_from_laycon(): - from flopy.utils.classic_to_mf6 import get_icelltype_from_laycon - # scalar confined (laycon=0) assert get_icelltype_from_laycon(0) == 0 # scalar convertible (laycon=1) assert get_icelltype_from_laycon(1) == 1 - # laycon=2 is head-dependent-T confined — maps to 0 + # laycon=2 is head-dependent-T confined assert get_icelltype_from_laycon(2) == 0 - # laycon=3 is fully convertible — maps to 1 + # laycon=3 is fully convertible assert get_icelltype_from_laycon(3) == 1 - # full 4-value array covering all cases + # array laycon = np.array([0, 1, 2, 3]) result = get_icelltype_from_laycon(laycon) assert np.array_equal(result, [0, 1, 0, 1]) - # shape is preserved + # make sure shape is preserved laycon_3d = np.array([[[0, 1], [2, 3]]]) result_3d = get_icelltype_from_laycon(laycon_3d) assert result_3d.shape == laycon_3d.shape assert np.array_equal(result_3d, [[[0, 1], [0, 1]]]) -def test_mfgrdfile_write_roundtrip(tmp_path): - from flopy.mf6.utils import MfGrdFile, get_structured_connectivity - - nlay, nrow, ncol = 2, 3, 4 - delr = np.ones(ncol) * 100.0 - delc = np.ones(nrow) * 50.0 - top = np.ones((nrow, ncol)) * 10.0 - botm = np.zeros((nlay, nrow, ncol)) - botm[0] = 5.0 - - ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) - - icelltype = np.zeros((nlay, nrow, ncol), dtype=np.int32) - icelltype[0] = 1 # top layer convertible - - grb_file = tmp_path / "test.dis.grb" - MfGrdFile.write_dis( - str(grb_file), - nlay, - nrow, - ncol, - delr, - delc, - top, - botm, - ia, - ja, - icelltype=icelltype, - ) - - grb = MfGrdFile(str(grb_file)) - - assert grb.nlay == nlay - assert grb.nrow == nrow - assert grb.ncol == ncol - assert grb.nodes == nlay * nrow * ncol - assert grb.nja == nja - - np.testing.assert_array_almost_equal(grb.delr, delr) - np.testing.assert_array_almost_equal(grb.delc, delc) - - # Build expected TOP for all cells (MF6 format) - # Layer 1: model top, Layer 2+: bottom of layer above - # In Fortran order (layer-interleaved): [L0[0,0], L1[0,0], L2[0,0], L0[0,1], ...] - top_expected = np.zeros(nlay * nrow * ncol) - top_flat = top.flatten(order="F") - botm_flat = botm.flatten(order="F") - - for i in range(nrow * ncol): - # Layer 0: use model top - top_expected[i * nlay] = top_flat[i] - # Layers 1+: use bottom of layer above - for k in range(1, nlay): - top_expected[i * nlay + k] = botm_flat[i * nlay + (k - 1)] - - np.testing.assert_array_almost_equal(grb.top, top_expected) - np.testing.assert_array_almost_equal(grb.bot, botm.flatten(order="F")) - - np.testing.assert_array_equal(grb.ia, ia) - np.testing.assert_array_equal(grb.ja, ja) - - icelltype_read = grb._datadict["ICELLTYPE"] - np.testing.assert_array_equal(icelltype_read, icelltype.flatten(order="F")) - - -def test_mfgrdfile_write_with_idomain(tmp_path): - from flopy.mf6.utils import MfGrdFile, get_structured_connectivity - - nlay, nrow, ncol = 1, 3, 3 - delr = np.ones(ncol) - delc = np.ones(nrow) - top = np.ones((nrow, ncol)) - botm = np.zeros((nlay, nrow, ncol)) - - # center cell inactive - idomain = np.ones((nlay, nrow, ncol), dtype=np.int32) - idomain[0, 1, 1] = 0 - - ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol, idomain) - - grb_file = tmp_path / "test_idomain.dis.grb" - MfGrdFile.write_dis( - str(grb_file), nlay, nrow, ncol, delr, delc, top, botm, ia, ja, idomain=idomain - ) - - grb = MfGrdFile(str(grb_file)) - idomain_read = grb._datadict["IDOMAIN"] - np.testing.assert_array_equal(idomain_read, idomain.flatten(order="F")) - - -def test_mfgrdfile_write_validation(): - from flopy.mf6.utils import MfGrdFile, get_structured_connectivity - - nlay, nrow, ncol = 1, 2, 2 - delr = np.ones(ncol) - delc = np.ones(nrow) - top = np.ones((nrow, ncol)) - botm = np.zeros((nlay, nrow, ncol)) - ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) - - with pytest.raises(ValueError, match="delr length"): - MfGrdFile.write_dis( - "test.grb", - nlay, - nrow, - ncol, - np.ones(3), # Wrong length - delc, - top, - botm, - ia, - ja, - ) - - with pytest.raises(ValueError, match="ia length"): - MfGrdFile.write_dis( - "test.grb", - nlay, - nrow, - ncol, - delr, - delc, - top, - botm, - np.ones(10), # Wrong length (should be ncells + 1 = 5) - ja, - ) - - -def test_nwt_to_mf6_converter_init(function_tmpdir): - import numpy as np - - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile - - nlay, nrow, ncol = 2, 3, 4 - delr = np.ones(ncol) * 100.0 - delc = np.ones(nrow) * 100.0 - top = np.ones((nrow, ncol)) * 10.0 - botm = np.zeros((nlay, nrow, ncol)) - botm[0] = 5.0 - laytyp = np.array([1, 0]) # Top layer convertible - - hds_file = function_tmpdir / "test.hds" - cbc_file = function_tmpdir / "test.cbc" - - head_data = np.ones((nlay, nrow, ncol)) * 8.0 - HeadFile.write( - str(hds_file), - {(1, 1): head_data}, - nlay=nlay, - nrow=nrow, - ncol=ncol, - ) - - from flopy.mf6.utils import get_structured_connectivity, get_structured_faceflows - - ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) - flowja = np.ones(nja) * 0.5 - - frf, fff, flf = get_structured_faceflows( - flowja, grb_file=None, ia=ia, ja=ja, nlay=nlay, nrow=nrow, ncol=ncol - ) - - bud_data = [ - { - "data": frf.flatten(order="F"), - "kstp": 1, - "kper": 1, - "text": "FLOW RIGHT FACE", - "imeth": 1, - }, - { - "data": fff.flatten(order="F"), - "kstp": 1, - "kper": 1, - "text": "FLOW FRONT FACE", - "imeth": 1, - }, - { - "data": flf.flatten(order="F"), - "kstp": 1, - "kper": 1, - "text": "FLOW LOWER FACE", - "imeth": 1, - }, - ] - - CellBudgetFile.write( - str(cbc_file), - bud_data, - nlay=nlay, - nrow=nrow, - ncol=ncol, - ) - - converter = ClassicMfToMf6Converter( - str(hds_file), - str(cbc_file), - nlay, - nrow, - ncol, - delr, - delc, - top, - botm, - laytyp, - model_ws=function_tmpdir, - ) - - assert converter.nlay == nlay - assert converter.nrow == nrow - assert converter.ncol == ncol - assert converter.ncells == nlay * nrow * ncol - assert len(converter.times) == 1 - assert len(converter.kstpkper) == 1 - - assert converter.icelltype.shape == (nlay,) - assert np.array_equal(converter.icelltype, [1, 0]) - assert converter.icelltype_3d.shape == (nlay, nrow, ncol) - assert np.all(converter.icelltype_3d[0] == 1) - assert np.all(converter.icelltype_3d[1] == 0) - - -def test_nwt_to_mf6_converter_convert(function_tmpdir): - import numpy as np - - from flopy.mf6.utils import ( - MfGrdFile, - get_structured_connectivity, - get_structured_faceflows, - ) - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile - +def test_converter_synthetic(function_tmpdir): nlay, nrow, ncol = 2, 3, 4 delr = np.ones(ncol) * 100.0 delc = np.ones(nrow) * 100.0 @@ -489,7 +158,7 @@ def test_nwt_to_mf6_converter_convert(function_tmpdir): @requires_exe("mfnwt") @pytest.mark.slow -def test_nwt_to_mf6_watertable_model(function_tmpdir): +def test_converter_mfnwt(function_tmpdir): from flopy.mf6.utils import MfGrdFile from flopy.modflow import ( Modflow, @@ -504,7 +173,7 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): from flopy.utils import ClassicMfToMf6Converter from flopy.utils.binaryfile import CellBudgetFile, HeadFile - modelname = "watertable" + modelname = "mfnwt" nlay, nrow, ncol = 1, 1, 100 delr = 50.0 @@ -569,7 +238,7 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): mf.write_input() success, _ = mf.run_model(silent=True) - assert success, "NWT model run failed" + assert success hds_file = function_tmpdir / f"{modelname}.hds" cbc_file = function_tmpdir / f"{modelname}.cbc" @@ -632,7 +301,7 @@ def test_nwt_to_mf6_watertable_model(function_tmpdir): @requires_exe("mfnwt") @pytest.mark.slow -def test_nwt_to_mf6_multilayer_model(function_tmpdir): +def test_converter_mfnwt_multilayer(function_tmpdir): from flopy.mf6.utils import MfGrdFile from flopy.modflow import ( Modflow, @@ -647,7 +316,7 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): from flopy.utils import ClassicMfToMf6Converter from flopy.utils.binaryfile import CellBudgetFile, HeadFile - modelname = "multilayer" + modelname = "mfnwt_multil" # 3 layers: top convertible, middle/bottom confined nlay, nrow, ncol = 3, 10, 10 @@ -797,11 +466,6 @@ def test_nwt_to_mf6_multilayer_model(function_tmpdir): assert len(sat_data) > 0, "Could not read DATA-SAT" -# --------------------------------------------------------------------------- -# Fixtures for example-data models -# --------------------------------------------------------------------------- - - @pytest.fixture def freyberg_multilayer_path(example_data_path): return example_data_path / "freyberg_multilayer_transient" @@ -812,29 +476,19 @@ def mf2005_freyberg_path(example_data_path): return example_data_path / "freyberg" -# --------------------------------------------------------------------------- -# Integration tests against real pre-computed model output -# --------------------------------------------------------------------------- - - @pytest.mark.slow -def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_tmpdir): +def test_converter_freyberg_multilayer(freyberg_multilayer_path, function_tmpdir): """ - Convert the pre-computed freyberg_multilayer_transient NWT outputs to MF6 + Convert the freyberg_multilayer_transient NWT outputs to MF6 and verify that the face-flow roundtrip is exact. - This test requires no executable — it uses the binary output files already - committed to the repo. It checks: + This test checks: - GRB file has the correct grid dimensions - Head values are preserved for a sample of time steps (first, middle, last) - - get_structured_faceflows(flowja_converted) recovers the original NWT + - get_structured_faceflows(flowja_converted) recovers the original FLOW RIGHT/FRONT/LOWER FACE values (interior faces only) - DATA-SAT is present for sampled time steps """ - from flopy.mf6.utils import MfGrdFile, get_structured_faceflows - from flopy.modflow import Modflow - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile model_ws = freyberg_multilayer_path @@ -866,7 +520,6 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t output_dir = function_tmpdir / "mf6" result = converter.convert(str(output_dir)) - # --- GRB dimensions --- grb = MfGrdFile(str(result["grb"])) assert grb.nlay == nlay assert grb.nrow == nrow @@ -879,12 +532,12 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t kstpkper_list = hds_orig.get_kstpkper() - # --- Verify a sample of time steps (first, middle, last) --- + # check sample of time steps (first, middle, last) check_indices = [0, len(kstpkper_list) // 2, len(kstpkper_list) - 1] for check_idx in check_indices: kstpkper = kstpkper_list[check_idx] - # Head values are preserved exactly + # check head values head_orig = hds_orig.get_data(kstpkper=kstpkper) head_mf6 = hds_mf6.get_data(kstpkper=kstpkper) np.testing.assert_array_almost_equal( @@ -894,7 +547,7 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t err_msg=f"Head mismatch at kstpkper={kstpkper}", ) - # Face-flow roundtrip: faceflows → flowja → faceflows should recover + # check face-flow roundtrip: faceflows → flowja → faceflows should recover # the original values on interior faces. frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] @@ -924,29 +577,21 @@ def test_classic_to_mf6_freyberg_multilayer(freyberg_multilayer_path, function_t err_msg=f"FLOW LOWER FACE mismatch at kstpkper={kstpkper}", ) - # DATA-SAT is present + # check DATA-SAT is present sat = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) assert len(sat) > 0, f"DATA-SAT missing at kstpkper={kstpkper}" @requires_exe("mf2005") @pytest.mark.slow -def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): +def test_converter_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): """ Run the single-layer steady-state Freyberg MODFLOW-2005 model, convert its outputs to MF6, and verify correctness. - Tests that the converter works with the LPF package (vs UPW for NWT), - confirming ClassicMfToMf6Converter handles both code paths. + Tests that the converter works with the LPF package (vs UPW for NWT). """ - import shutil - - from flopy.mf6.utils import MfGrdFile, get_structured_faceflows - from flopy.modflow import Modflow - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile - # Copy model to a writable temp directory and run it run_ws = function_tmpdir / "mf2005" shutil.copytree(mf2005_freyberg_path, run_ws) @@ -985,14 +630,14 @@ def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): kstpkper = hds_orig.get_kstpkper()[0] - # Head roundtrip + # check head roundtrip np.testing.assert_array_almost_equal( hds_mf6.get_data(kstpkper=kstpkper), hds_orig.get_data(kstpkper=kstpkper), decimal=5, ) - # Face-flow roundtrip (single layer, so no lower face) + # check face-flow roundtrip (single layer, so no lower face) frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] @@ -1009,16 +654,11 @@ def test_classic_to_mf6_freyberg_mf2005(mf2005_freyberg_path, function_tmpdir): @requires_exe("mf2000") @pytest.mark.slow -def test_classic_to_mf6_mf2000_watertable(function_tmpdir): +def test_converter_mf2000_1layer(function_tmpdir): """ Build and run a 1-layer unconfined watertable MODFLOW-2000 model, convert its outputs to MF6 format, and verify head, face-flow, and saturation. - - MODFLOW-2000 produces compact budget files with the same binary format as - MODFLOW-NWT and MODFLOW-2005. This test confirms ClassicMfToMf6Converter - handles MF-2000 output identically to later variants. """ - from flopy.mf6.utils import MfGrdFile, get_structured_faceflows from flopy.modflow import ( Modflow, ModflowBas, @@ -1029,10 +669,8 @@ def test_classic_to_mf6_mf2000_watertable(function_tmpdir): ModflowPcg, ModflowRch, ) - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile - modelname = "mf2000_watertable" + modelname = "mf2000_1layer" nlay, nrow, ncol = 1, 1, 100 delr, delc = 50.0, 1.0 h1, h2 = 20.0, 11.0 @@ -1109,13 +747,13 @@ def test_classic_to_mf6_mf2000_watertable(function_tmpdir): kstpkper = hds_orig.get_kstpkper()[0] - # Head roundtrip + # check head roundtrip head_orig = hds_orig.get_data(kstpkper=kstpkper) np.testing.assert_array_almost_equal( hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 ) - # Face-flow roundtrip (1 row, 1 layer → only right face flows exist) + # check face-flow roundtrip (1 row, 1 layer → only right face flows exist) frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] frf_rt, _, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) @@ -1123,7 +761,7 @@ def test_classic_to_mf6_mf2000_watertable(function_tmpdir): frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 ) - # Saturation: sat = clamp(head / (top - botm), 0, 1) + # check saturation: sat = clamp(head / (top - botm), 0, 1) sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) assert len(sat_data) > 0, "DATA-SAT missing" expected_sat = np.clip(head_orig.flatten() / (top - botm), 0.0, 1.0) @@ -1139,7 +777,7 @@ def test_classic_to_mf6_mf2000_watertable(function_tmpdir): "mf6/test/test001h_rch_array3", # 1 layer, 4 stress periods ], ) -def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): +def test_converter_mf2005_testmodels(model_key, function_tmpdir): """ Fetch a MODFLOW-2005/LPF companion model from the devtools registry, run it, convert outputs to MF6, and verify head and face-flow roundtrip. @@ -1148,27 +786,22 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): test case. Running these models exercises the LPF code path across a variety of boundary conditions and confirms multi-stress-period conversion. """ - import shutil + from modflow_devtools import models - from modflow_devtools.models import copy_to - - from flopy.mf6.utils import get_structured_faceflows - from flopy.modflow import Modflow - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile + from flopy.modflow import ModflowOc - # Download the parent MF6 model; the mf2005/ companion is a subdirectory. + # fetch the model model_ws = function_tmpdir / "model" - copy_to(str(model_ws), model_key) + models.copy_to(str(model_ws), model_key) mf2005_ws = model_ws / "mf2005" assert mf2005_ws.exists(), f"mf2005/ subdir not found for {model_key}" - # Locate the namefile (there should be exactly one .nam in mf2005/) + # find namefile nam_files = list(mf2005_ws.glob("*.nam")) assert len(nam_files) == 1, f"Expected 1 .nam file in mf2005/, got {nam_files}" nam_name = nam_files[0].name - # Run mf2005 + # load model mf = Modflow.load( nam_name, model_ws=str(mf2005_ws), @@ -1176,8 +809,7 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): check=False, ) - # The companion models may not have "save budget" in their OC. Force it so - # that the CBC file is populated for the converter. + # make sure output saving is enabled nper = int(mf.dis.nper) nstp = mf.dis.nstp.array sp_data = { @@ -1185,7 +817,6 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): for kper in range(nper) for kstp in range(nstp[kper]) } - from flopy.modflow import ModflowOc ModflowOc(mf, stress_period_data=sp_data, compact=True) mf.write_input() @@ -1193,7 +824,6 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): success, _ = mf.run_model(silent=True) assert success, f"mf2005 run failed for {model_key}" - # Locate binary outputs written by the OC package hds_path = next(mf2005_ws.glob("*.hds"), None) cbc_path = next(mf2005_ws.glob("*.cbc"), None) assert hds_path is not None, "No .hds file produced" @@ -1227,7 +857,7 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): bud_mf6 = CellBudgetFile(str(result["bud"])) for kstpkper in hds_orig.get_kstpkper(): - # Head roundtrip + # check head roundtrip np.testing.assert_array_almost_equal( hds_mf6.get_data(kstpkper=kstpkper), hds_orig.get_data(kstpkper=kstpkper), @@ -1235,7 +865,7 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): err_msg=f"Head mismatch at kstpkper={kstpkper}", ) - # Face-flow roundtrip (check only terms that the model actually wrote) + # check face-flow roundtrip (only terms the model actually wrote) cbc_texts = [ t.decode().strip() if isinstance(t, bytes) else t.strip() for t in cbc_orig.textlist @@ -1276,15 +906,11 @@ def test_classic_to_mf6_mf2005_testmodels(model_key, function_tmpdir): @requires_exe("mf2005") @pytest.mark.slow -def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): +def test_converter_mf2005_multilayer(function_tmpdir): """ Build and run a 3-layer MODFLOW-2005/LPF model (1 convertible + 2 confined), convert its outputs to MF6, and verify head, face-flow, and saturation. - - Complements test_classic_to_mf6_freyberg_mf2005 (single layer) and the - NWT multilayer tests by exercising the LPF code path with nlay > 1. """ - from flopy.mf6.utils import MfGrdFile, get_structured_faceflows from flopy.modflow import ( Modflow, ModflowBas, @@ -1295,8 +921,6 @@ def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): ModflowRch, ModflowWel, ) - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile modelname = "mf2005_multilayer" nlay, nrow, ncol = 3, 10, 10 @@ -1353,13 +977,11 @@ def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): converter = ClassicMfToMf6Converter.from_model(mf, hds_path, cbc_path) result = converter.convert(str(function_tmpdir / "mf6_output")) - # GRB dimensions grb = MfGrdFile(str(result["grb"])) assert grb.nlay == nlay assert grb.nrow == nrow assert grb.ncol == ncol - # ICELLTYPE pattern (Fortran interleaved: layer varies fastest) icelltype = grb._datadict["ICELLTYPE"] assert np.all(icelltype[0::nlay] == 1), "Layer 1 should be convertible" assert np.all(icelltype[1::nlay] == 0), "Layer 2 should be confined" @@ -1373,12 +995,12 @@ def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): kstpkper = hds_orig.get_kstpkper()[0] head_orig = hds_orig.get_data(kstpkper=kstpkper) - # Head roundtrip + # check head roundtrip np.testing.assert_array_almost_equal( hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 ) - # Face-flow roundtrip (3 layers → include lower face) + # head face-flow roundtrip (3 layers → include lower face) frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] fff_orig = cbc_orig.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] flf_orig = cbc_orig.get_data(text="FLOW LOWER FACE", kstpkper=kstpkper)[0] @@ -1396,18 +1018,20 @@ def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): flf_rt[:-1, :, :], flf_orig[:-1, :, :], decimal=5 ) - # Saturation: confined layers = 1.0; convertible layer = (head - botm) / thickness + # check saturation: + # confined layers = 1.0; + # convertible layer = (head - botm) / thickness sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) assert len(sat_data) > 0, "DATA-SAT missing" - sat_flat = sat_data[0]["q"] # Fortran order: layer varies fastest - # Confined layers (1 and 2, 0-indexed) should all be 1.0 + sat_flat = sat_data[0]["q"] + # confined layers (1 and 2, 0-indexed) should all be 1.0 np.testing.assert_array_equal( sat_flat[1::nlay], 1.0, err_msg="Layer 2 (confined) sat != 1" ) np.testing.assert_array_equal( sat_flat[2::nlay], 1.0, err_msg="Layer 3 (confined) sat != 1" ) - # Convertible layer sat should match (head - botm[0]) / (top - botm[0]) + # convertible layer sat should match (head - botm[0]) / (top - botm[0]) head_layer0 = head_orig[0].flatten(order="F") top_layer0 = top.flatten(order="F") bot_layer0 = botm[0].flatten(order="F") @@ -1421,19 +1045,13 @@ def test_classic_to_mf6_mf2005_multilayer(function_tmpdir): @requires_exe("mf2000") @pytest.mark.slow -def test_classic_to_mf6_mf2000_bcf(function_tmpdir): +def test_converter_mf2000_bcf(function_tmpdir): """ Build and run a MODFLOW-2000 model with the BCF package (laycon=3, - fully convertible), convert via from_model(), and verify that - get_icelltype_from_laycon() correctly maps laycon → ICELLTYPE and that - head, face-flow, and saturation all round-trip correctly. - - This test exercises the BCF code path in from_model(), which differs - from the LPF/UPW path: laycon values (0/1/2/3) are first converted to - ICELLTYPE (0/1) via get_icelltype_from_laycon(), and hdry is read from - the BCF package (default -1e30, not -999.0). + fully convertible), then perform the conversion with from_model(), + and verify that get_icelltype_from_laycon() correctly maps laycon + to cell type and head/face-flows/saturation round-trip correctly. """ - from flopy.mf6.utils import get_structured_faceflows from flopy.modflow import ( Modflow, ModflowBas, @@ -1444,8 +1062,6 @@ def test_classic_to_mf6_mf2000_bcf(function_tmpdir): ModflowPcg, ModflowRch, ) - from flopy.utils import ClassicMfToMf6Converter - from flopy.utils.binaryfile import CellBudgetFile, HeadFile modelname = "mf2000_bcf" nlay, nrow, ncol = 1, 1, 100 @@ -1454,7 +1070,6 @@ def test_classic_to_mf6_mf2000_bcf(function_tmpdir): top, botm = 25.0, 0.0 hk = 50.0 - # Linear initial heads to avoid convergence failures with PCG and unconfined cells strt = np.zeros((nlay, nrow, ncol)) strt[0, 0, :] = np.linspace(h1, h2, ncol) @@ -1528,12 +1143,12 @@ def test_classic_to_mf6_mf2000_bcf(function_tmpdir): kstpkper = hds_orig.get_kstpkper()[0] head_orig = hds_orig.get_data(kstpkper=kstpkper) - # Head roundtrip + # check head roundtrip np.testing.assert_array_almost_equal( hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 ) - # Face-flow roundtrip (1 row, 1 layer → only right face flows exist) + # check face-flow roundtrip (1 row, 1 layer → only right face flows exist) frf_orig = cbc_orig.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] flowja = bud_mf6.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)[0] frf_rt, _, _ = get_structured_faceflows(flowja, grb_file=str(result["grb"])) @@ -1541,7 +1156,7 @@ def test_classic_to_mf6_mf2000_bcf(function_tmpdir): frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 ) - # Saturation (laycon=3 → convertible → sat = clamp(head / 25, 0, 1)) + # check saturation (laycon=3 → convertible → sat = clamp(head / 25, 0, 1)) sat_data = bud_mf6.get_data(text="DATA-SAT", kstpkper=kstpkper) assert len(sat_data) > 0, "DATA-SAT missing" active = head_orig.flatten() > -1e29 # exclude hdry/hnoflo From ac7c6ffa1f51ea8d242bba83a24b684827497179 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 15:38:06 -0700 Subject: [PATCH 09/10] more cleanup --- flopy/utils/classic_to_mf6.py | 186 +++++++++++++--------------------- 1 file changed, 69 insertions(+), 117 deletions(-) diff --git a/flopy/utils/classic_to_mf6.py b/flopy/utils/classic_to_mf6.py index bae8f29f30..6195d5a588 100644 --- a/flopy/utils/classic_to_mf6.py +++ b/flopy/utils/classic_to_mf6.py @@ -166,7 +166,7 @@ def get_icelltype_from_laycon(laycon): class ClassicMfToMf6Converter: """ - Convert classic MODFLOW binary outputs to MODFLOW 6 binary format. + Convert classic MODFLOW compact binary files to MODFLOW 6 binary format. Reads head and cell-budget files produced by any structured-grid MODFLOW variant with compact budget output (FLOW RIGHT FACE / FLOW FRONT FACE / @@ -174,13 +174,12 @@ class ClassicMfToMf6Converter: file (``FLOW-JA-FACE``, ``DATA-SAT``), and binary grid record (GRB) required by PRT and other MF6 post-processors. - Confirmed compatible variants: MODFLOW-NWT (UPW), MODFLOW-2005 (LPF), - MODFLOW-2000 (LPF or BCF). + Compatible with: MODFLOW-NWT (UPW), MODFLOW-2005 (LPF), MODFLOW-2000 (LPF or BCF). The easiest way to construct a converter from a loaded model is :meth:`from_model`, which autodetects the flow package and extracts all - required parameters automatically. Construct directly only when a model - object is not available (e.g. you only have the binary files on disk). + required parameters automatically. A converter can also be constructed + directly if a model is not available (e.g. you have binary output files). Parameters ---------- @@ -469,7 +468,6 @@ def convert( - 'bud': Path to budget file """ from ..mf6.utils import MfGrdFile - from .binaryfile import CellBudgetFile, HeadFile # Create output directory output_dir = Path(output_dir) @@ -523,22 +521,20 @@ def convert( return {"grb": grb_path, "hds": hds_path, "bud": bud_path} - def _build_header_lookup(self): + def _build_budget_header_lookup(self): """ - Build a dict mapping (kstp, kper) -> header record from the head file. - - Returns - ------- - dict - Keys are (kstp, kper) tuples, values are header records with - fields including 'pertim' and 'totim'. + Build a dict mapping 0-based (kstp, kper) -> first CBC header record + for that time step. The compact budget header stores delt, pertim, + and totim directly, so no arithmetic is needed to recover them. """ - # get_kstpkper() returns 0-based indices; recordarray stores 1-based - # file values. Subtract 1 so the lookup keys match self.kstpkper. - return { - (int(rec["kstp"]) - 1, int(rec["kper"]) - 1): rec - for rec in self.hds_obj.recordarray - } + seen = set() + result = {} + for rec in self.cbc_obj.recordarray: + key = (int(rec["kstp"]) - 1, int(rec["kper"]) - 1) + if key not in seen: + result[key] = rec + seen.add(key) + return result def _write_budget(self, filename, precision, verbose): """Write budget file with FLOW-JA-FACE, DATA-SPDIS, and DATA-SAT.""" @@ -549,114 +545,83 @@ def _write_budget(self, filename, precision, verbose): if verbose: print(f" Processing {len(self.kstpkper)} time steps...") - # We'll write three terms per time step: - # 1. FLOW-JA-FACE (imeth=1, array) - # 2. DATA-SPDIS (imeth=6, list with qx, qy, qz) - # 3. DATA-SAT (imeth=6, list) - - # Build list of records records = [] + header_lookup = self._build_budget_header_lookup() - header_lookup = self._build_header_lookup() - - for idx, kstpkper in enumerate(self.kstpkper): + for kstpkper in self.kstpkper: kstp, kper = kstpkper rec = header_lookup[kstpkper] totim = float(rec["totim"]) pertim = float(rec["pertim"]) - - # delt: for the first time step within a period pertim equals delt; - # for subsequent steps subtract the previous step's pertim. - # kstp is 0-based here (from get_kstpkper()). - if kstp == 0: - delt = pertim - else: - prev_rec = header_lookup.get((kstp - 1, kper)) - delt = pertim - float(prev_rec["pertim"]) if prev_rec else pertim + delt = float(rec["delt"]) if verbose: print(f" Processing time step {kstpkper}...") - # Get head head = self.hds_obj.get_data(kstpkper=kstpkper) - # Get face flows if verbose: print( f" Available budget terms: " f"{self.cbc_obj.get_unique_record_names()}" ) - # Check which face flows are available - # For 1D/2D models, not all face flows may exist + # for 1D/2D models, some directions may be missing available_terms = [ t.decode().strip() for t in self.cbc_obj.get_unique_record_names() ] - try: - # FLOW RIGHT FACE (required for X-direction flow) - if "FLOW RIGHT FACE" in available_terms: - frf_data = self.cbc_obj.get_data( - text="FLOW RIGHT FACE", kstpkper=kstpkper - ) - if verbose: - print( - f" FLOW RIGHT FACE: {type(frf_data)}, " - f"len={len(frf_data) if frf_data else 0}" - ) - frf = frf_data[0] if frf_data and len(frf_data) > 0 else None - else: - frf = None - - # FLOW FRONT FACE (required for Y-direction flow) - if "FLOW FRONT FACE" in available_terms: - fff_data = self.cbc_obj.get_data( - text="FLOW FRONT FACE", kstpkper=kstpkper + if "FLOW RIGHT FACE" in available_terms: + frf_data = self.cbc_obj.get_data( + text="FLOW RIGHT FACE", kstpkper=kstpkper + ) + if verbose: + print( + f" FLOW RIGHT FACE: {type(frf_data)}, " + f"len={len(frf_data) if frf_data else 0}" ) - if verbose: - print( - f" FLOW FRONT FACE: {type(fff_data)}, " - f"len={len(fff_data) if fff_data else 0}" - ) - fff = fff_data[0] if fff_data and len(fff_data) > 0 else None - else: - fff = None - - # FLOW LOWER FACE (required for Z-direction flow) - if "FLOW LOWER FACE" in available_terms: - flf_data = self.cbc_obj.get_data( - text="FLOW LOWER FACE", kstpkper=kstpkper + frf = frf_data[0] if frf_data and len(frf_data) > 0 else None + else: + frf = None + + if "FLOW FRONT FACE" in available_terms: + fff_data = self.cbc_obj.get_data( + text="FLOW FRONT FACE", kstpkper=kstpkper + ) + if verbose: + print( + f" FLOW FRONT FACE: {type(fff_data)}, " + f"len={len(fff_data) if fff_data else 0}" ) - if verbose: - print( - f" FLOW LOWER FACE: {type(flf_data)}, " - f"len={len(flf_data) if flf_data else 0}" - ) - flf = flf_data[0] if flf_data and len(flf_data) > 0 else None - else: - flf = None - - # Validate at least one face flow exists - if frf is None and fff is None and flf is None: - raise ValueError("No face flows found in budget file") - - # Create zero arrays for missing face flows - # For 1D/2D models, not all directions have flow - shape_3d = (self.nlay, self.nrow, self.ncol) - if frf is None: - frf = np.zeros(shape_3d, dtype=np.float64) - if fff is None: - fff = np.zeros(shape_3d, dtype=np.float64) - if flf is None: - flf = np.zeros(shape_3d, dtype=np.float64) - - except Exception as e: + fff = fff_data[0] if fff_data and len(fff_data) > 0 else None + else: + fff = None + + if "FLOW LOWER FACE" in available_terms: + flf_data = self.cbc_obj.get_data( + text="FLOW LOWER FACE", kstpkper=kstpkper + ) if verbose: - print(f" Warning: Could not read face flows: {e}") - print(f" Skipping time step {kstpkper}") - continue + print( + f" FLOW LOWER FACE: {type(flf_data)}, " + f"len={len(flf_data) if flf_data else 0}" + ) + flf = flf_data[0] if flf_data and len(flf_data) > 0 else None + else: + flf = None + + if frf is None and fff is None and flf is None: + raise ValueError("No face flows found in budget file") + + # create zero arrays for missing face flows + shape_3d = (self.nlay, self.nrow, self.ncol) + if frf is None: + frf = np.zeros(shape_3d, dtype=np.float64) + if fff is None: + fff = np.zeros(shape_3d, dtype=np.float64) + if flf is None: + flf = np.zeros(shape_3d, dtype=np.float64) - # 1. Convert to FLOW-JA-FACE flowja = get_structured_flowja( (frf, fff, flf), ia=self.ia, @@ -665,7 +630,6 @@ def _write_budget(self, filename, precision, verbose): nrow=self.nrow, ncol=self.ncol, ) - records.append( { "data": flowja, @@ -679,25 +643,19 @@ def _write_budget(self, filename, precision, verbose): } ) - # 2. DATA-SPDIS (specific discharge) - SKIPPED for now # TODO: get_specific_discharge() requires a model object and cannot - # easily be driven from binary files alone. PRT can reconstruct - # specific discharge from FLOW-JA-FACE if needed. + # easily be driven from binary files alone. PRT doesn't need SPDIS, + # but other models might, so skip it for now and come back to this. - # 3. Calculate saturation sat = get_saturation( head, self.top, self.botm, self.icelltype_3d, self.hdry, self.hnoflo ) - # Build list data for DATA-SAT sat_flat = sat.flatten(order="F") active_sat = ~np.isnan(sat_flat) nlist_sat = np.sum(active_sat) - if nlist_sat > 0: nodes_sat = np.arange(self.ncells)[active_sat] + 1 # 1-based - - # Create structured array for imeth=6 dtype = np.dtype( [ ("node", np.int32), @@ -709,7 +667,6 @@ def _write_budget(self, filename, precision, verbose): sat_data["node"] = nodes_sat sat_data["node2"] = nodes_sat sat_data["sat"] = sat_flat[active_sat] - records.append( { "data": sat_data, @@ -724,7 +681,6 @@ def _write_budget(self, filename, precision, verbose): } ) - # Write all records if verbose: print(f" Writing {len(records)} budget records...") @@ -747,7 +703,3 @@ def __repr__(self): f" time_steps={len(self.times)}\n" f")" ) - - -#: Backward-compatible alias. New code should use ClassicMfToMf6Converter. -NwtToMf6Converter = ClassicMfToMf6Converter From 19cd49c0c186b03762fdaea01eff60c688c036eb Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 16 Apr 2026 15:40:44 -0700 Subject: [PATCH 10/10] use setdefault, hoist available_terms outside loop, remove redundant conditions --- flopy/utils/__init__.py | 1 - flopy/utils/classic_to_mf6.py | 77 +++++++++++------------------------ 2 files changed, 23 insertions(+), 55 deletions(-) diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 3026f43bb4..dca51c0f3f 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -31,7 +31,6 @@ get_modflow = get_modflow_module.run_main from .classic_to_mf6 import ( ClassicMfToMf6Converter, - NwtToMf6Converter, # backward-compatible alias get_icelltype_from_laycon, get_icelltype_from_laytyp, ) diff --git a/flopy/utils/classic_to_mf6.py b/flopy/utils/classic_to_mf6.py index 6195d5a588..df651baab8 100644 --- a/flopy/utils/classic_to_mf6.py +++ b/flopy/utils/classic_to_mf6.py @@ -527,13 +527,9 @@ def _build_budget_header_lookup(self): for that time step. The compact budget header stores delt, pertim, and totim directly, so no arithmetic is needed to recover them. """ - seen = set() result = {} for rec in self.cbc_obj.recordarray: - key = (int(rec["kstp"]) - 1, int(rec["kper"]) - 1) - if key not in seen: - result[key] = rec - seen.add(key) + result.setdefault((int(rec["kstp"]) - 1, int(rec["kper"]) - 1), rec) return result def _write_budget(self, filename, precision, verbose): @@ -548,6 +544,13 @@ def _write_budget(self, filename, precision, verbose): records = [] header_lookup = self._build_budget_header_lookup() + # for 1D/2D models, some face-flow directions may be absent + available_terms = { + t.decode().strip() for t in self.cbc_obj.get_unique_record_names() + } + if verbose: + print(f" Available budget terms: {available_terms}") + for kstpkper in self.kstpkper: kstp, kper = kstpkper rec = header_lookup[kstpkper] @@ -560,55 +563,21 @@ def _write_budget(self, filename, precision, verbose): head = self.hds_obj.get_data(kstpkper=kstpkper) - if verbose: - print( - f" Available budget terms: " - f"{self.cbc_obj.get_unique_record_names()}" - ) - - # for 1D/2D models, some directions may be missing - available_terms = [ - t.decode().strip() for t in self.cbc_obj.get_unique_record_names() - ] - - if "FLOW RIGHT FACE" in available_terms: - frf_data = self.cbc_obj.get_data( - text="FLOW RIGHT FACE", kstpkper=kstpkper - ) - if verbose: - print( - f" FLOW RIGHT FACE: {type(frf_data)}, " - f"len={len(frf_data) if frf_data else 0}" - ) - frf = frf_data[0] if frf_data and len(frf_data) > 0 else None - else: - frf = None - - if "FLOW FRONT FACE" in available_terms: - fff_data = self.cbc_obj.get_data( - text="FLOW FRONT FACE", kstpkper=kstpkper - ) - if verbose: - print( - f" FLOW FRONT FACE: {type(fff_data)}, " - f"len={len(fff_data) if fff_data else 0}" - ) - fff = fff_data[0] if fff_data and len(fff_data) > 0 else None - else: - fff = None - - if "FLOW LOWER FACE" in available_terms: - flf_data = self.cbc_obj.get_data( - text="FLOW LOWER FACE", kstpkper=kstpkper - ) - if verbose: - print( - f" FLOW LOWER FACE: {type(flf_data)}, " - f"len={len(flf_data) if flf_data else 0}" - ) - flf = flf_data[0] if flf_data and len(flf_data) > 0 else None - else: - flf = None + frf = ( + self.cbc_obj.get_data(text="FLOW RIGHT FACE", kstpkper=kstpkper)[0] + if "FLOW RIGHT FACE" in available_terms + else None + ) + fff = ( + self.cbc_obj.get_data(text="FLOW FRONT FACE", kstpkper=kstpkper)[0] + if "FLOW FRONT FACE" in available_terms + else None + ) + flf = ( + self.cbc_obj.get_data(text="FLOW LOWER FACE", kstpkper=kstpkper)[0] + if "FLOW LOWER FACE" in available_terms + else None + ) if frf is None and fff is None and flf is None: raise ValueError("No face flows found in budget file")