diff --git a/autotest/test_nwt_to_mf6.py b/autotest/test_nwt_to_mf6.py new file mode 100644 index 000000000..0bde892f9 --- /dev/null +++ b/autotest/test_nwt_to_mf6.py @@ -0,0 +1,884 @@ +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 + + # 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) + # 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): + 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 + + 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 = ClassicMfToMf6Converter( + 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 ClassicMfToMf6Converter + 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 = 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, + ) + + 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 ClassicMfToMf6Converter + 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 = 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, + ) + + 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 + + +# --------------------------------------------------------------------------- +# 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/autotest/test_postprocessing.py b/autotest/test_postprocessing.py index d3fbbd2f6..41cb8c5f5 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,291 @@ 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): 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): no right or front neighbor at boundary + # connections: diagonal + left + back = 3 (full bidirectional connectivity) + node_corner = 0 * nrow * ncol + 2 * ncol + 2 + 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(): + 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}" + + # 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]] + 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(): + 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 + ) + + # 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] + + # 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(): + 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] + # 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(): + 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 45de130dd..3c527fd78 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 09c256454..99f9402db 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 4f219736d..2c6cbd8ad 100644 --- a/flopy/mf6/utils/postprocessing.py +++ b/flopy/mf6/utils/postprocessing.py @@ -3,6 +3,109 @@ 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 + ----- + 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. + + Examples + -------- + >>> import numpy as np + >>> from flopy.mf6.utils import get_structured_connectivity + >>> nlay, nrow, ncol = 2, 3, 3 + >>> ia, ja, nja = get_structured_connectivity(nlay, nrow, ncol) + >>> print(f"Total cells: {nlay * nrow * ncol}, connections: {nja}") + Total cells: 18, connections: 84 + """ + 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})" + ) + + # 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): + if idomain[k, i, j] <= 0: + 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 + + 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 + + def get_structured_faceflows( flowja, grb_file=None, @@ -52,19 +155,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 +172,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 +230,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 +243,168 @@ 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): + k, i, j = np.unravel_index(n, (nlay, nrow, ncol)) + if idomain[k, i, j] <= 0: + continue + + for ipos in range(ia[n], ia[n + 1]): + m = ja[ipos] + + if m == n: + flowja[ipos] = 0.0 # diagonal (residual; zero here) + continue + + km, im, jm = np.unravel_index(m, (nlay, nrow, ncol)) + + # 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] + + # 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 (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] + + # 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] + + # 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 + + +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 +447,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 +467,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 bd15a59b8..08e94145f 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -40,6 +40,11 @@ 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/nwt_to_mf6.py new file mode 100644 index 000000000..4d73a35d3 --- /dev/null +++ b/flopy/utils/nwt_to_mf6.py @@ -0,0 +1,564 @@ +""" +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 (UPW) and MODFLOW-2005 (LPF), LAYTYP indicates layer + properties: + - LAYTYP = 0: Confined, transmissivity constant + - 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 + + Examples + -------- + >>> import numpy as np + >>> from flopy.utils.nwt_to_mf6 import get_icelltype_from_laytyp + >>> laytyp = np.array([1, 0, -1]) # Top and bottom layers convertible + >>> icelltype = get_icelltype_from_laytyp(laytyp) + >>> print(icelltype) + [1 0 1] + """ + laytyp = np.atleast_1d(laytyp).astype(np.int32) + icelltype = np.where(laytyp != 0, 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. + + 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 head file (.hds) + cbc_file : str + Path to 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 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. + 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 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 + ... ) + >>> 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 _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 + + 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 .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 = [] + + header_lookup = self._build_header_lookup() + + for idx, kstpkper in enumerate(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 + + 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": int(kstp) + 1, + "kper": int(kper) + 1, + "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": int(kstp) + 1, + "kper": int(kper) + 1, + "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"{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 diff --git a/flopy/utils/postprocessing.py b/flopy/utils/postprocessing.py index ec9416626..d393ef034 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