diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index e3308ef9c..77e75b595 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 5d56e0320..314d5bce5 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 new file mode 100644 index 000000000..519f19879 --- /dev/null +++ b/autotest/test_classic_to_mf6.py @@ -0,0 +1,1164 @@ +import shutil + +import numpy as np +import pytest +from modflow_devtools.markers import requires_exe + +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(): + + # 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) 1 + assert np.array_equal(icelltype, [0, 1, 1, 1, 1]) + + +def test_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 + assert get_icelltype_from_laycon(2) == 0 + # laycon=3 is fully convertible + assert get_icelltype_from_laycon(3) == 1 + + # array + laycon = np.array([0, 1, 2, 3]) + result = get_icelltype_from_laycon(laycon) + assert np.array_equal(result, [0, 1, 0, 1]) + + # 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_converter_synthetic(function_tmpdir): + 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_converter_mfnwt(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 = "mfnwt" + + 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 + + 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.from_model(mf, hds_file, cbc_file) + + 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 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") +@pytest.mark.slow +def test_converter_mfnwt_multilayer(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 = "mfnwt_multil" + + # 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.from_model(mf, hds_file, cbc_file) + + 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 + 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" + + +@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" + + +@pytest.mark.slow +def test_converter_freyberg_multilayer(freyberg_multilayer_path, function_tmpdir): + """ + Convert the freyberg_multilayer_transient NWT outputs to MF6 + and verify that the face-flow roundtrip is exact. + + 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 + FLOW RIGHT/FRONT/LOWER FACE values (interior faces only) + - DATA-SAT is present for sampled time steps + """ + + 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 = 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() + + # 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] + + # 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( + head_mf6, + head_orig, + decimal=5, + err_msg=f"Head mismatch at kstpkper={kstpkper}", + ) + + # 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] + 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}", + ) + + # 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_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). + """ + + 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.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] + + # check head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), + hds_orig.get_data(kstpkper=kstpkper), + decimal=5, + ) + + # 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] + + 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_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. + """ + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowGhb, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowRch, + ) + + modelname = "mf2000_1layer" + 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] + + # 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 + ) + + # 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"])) + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 + ) + + # 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) + 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_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. + + 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. + """ + from modflow_devtools import models + + from flopy.modflow import ModflowOc + + # fetch the model + model_ws = function_tmpdir / "model" + 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}" + + # 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 + + # load model + mf = Modflow.load( + nam_name, + model_ws=str(mf2005_ws), + exe_name="mf2005", + check=False, + ) + + # make sure output saving is enabled + 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]) + } + + 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}" + + 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), + 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, + ) + 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"])) + + for kstpkper in hds_orig.get_kstpkper(): + # check 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}", + ) + + # 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 + ] + 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_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. + """ + from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowRch, + ModflowWel, + ) + + 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 = MfGrdFile(str(result["grb"])) + assert grb.nlay == nlay + assert grb.nrow == nrow + assert grb.ncol == ncol + + 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) + + # check head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 + ) + + # 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] + 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 + ) + 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 + ) + + # 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"] + # 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_converter_mf2000_bcf(function_tmpdir): + """ + Build and run a MODFLOW-2000 model with the BCF package (laycon=3, + 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.modflow import ( + Modflow, + ModflowBas, + ModflowBcf, + ModflowDis, + ModflowGhb, + ModflowOc, + ModflowPcg, + ModflowRch, + ) + + 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 + + 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) + + # check head roundtrip + np.testing.assert_array_almost_equal( + hds_mf6.get_data(kstpkper=kstpkper), head_orig, decimal=5 + ) + + # 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"])) + np.testing.assert_array_almost_equal( + frf_rt[:, :, :-1], frf_orig[:, :, :-1], decimal=5 + ) + + # 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 + 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/autotest/test_postprocessing.py b/autotest/test_postprocessing.py index d3fbbd2f6..dfaa98cb4 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,293 @@ 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..dca51c0f3 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -29,6 +29,11 @@ from .formattedfile import FormattedHeadFile get_modflow = get_modflow_module.run_main +from .classic_to_mf6 import ( + ClassicMfToMf6Converter, + get_icelltype_from_laycon, + get_icelltype_from_laytyp, +) from .gridintersect import GridIntersect from .mflistfile import ( Mf6ListBudget, diff --git a/flopy/utils/classic_to_mf6.py b/flopy/utils/classic_to_mf6.py new file mode 100644 index 000000000..df651baab --- /dev/null +++ b/flopy/utils/classic_to_mf6.py @@ -0,0 +1,674 @@ +""" +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, 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: compact MODFLOW vs MODFLOW 6 +-------------------------------------------------------- +**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)``. + +**Budget files (.cbc / .bud)** + 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``); 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-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 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 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 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 + flow models create a binary grid file (``.dis.grb``) that encodes + the IA/JA connectivity, cell geometry, convertibility, and IDOMAIN. +""" + +import shutil +from pathlib import Path + +import numpy as np + + +def get_icelltype_from_laytyp(laytyp): + """ + Convert classic MODFLOW LAYTYP values to MF6 ICELLTYPE values. + + Parameters + ---------- + laytyp : array_like + 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) + + Returns + ------- + icelltype : ndarray + Cell type array for MF6 (same shape as input): + - 0: Confined + - 1: Convertible + + Notes + ----- + In MODFLOW-NWT (UPW) and MODFLOW-2005/2000 (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 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 + + +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 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 / + 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. + + 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. A converter can also be constructed + directly if a model is not available (e.g. you have binary output files). + + Parameters + ---------- + 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. + nrow : int + Number of rows. + ncol : int + 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 + 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 + 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 + 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 + >>> laytyp = np.array([1, 0, 0]) # top layer convertible + >>> converter = ClassicMfToMf6Converter( + ... 'model.hds', 'model.cbc', + ... nlay, nrow, ncol, + ... 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') + """ + + 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() + + @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, + grb_name="gwf.grb", + hds_name="gwf.hds", + bud_name="gwf.bud", + precision="double", + verbose=False, + ): + """ + Convert classic MODFLOW binary outputs to MF6 format. + + Parameters + ---------- + output_dir : str or PathLike + 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 + 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 + + # 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 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}") + 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, + ) + + # Copy head file verbatim — the on-disk format is identical between + # classic MODFLOW and MF6, so no conversion is needed. + if verbose: + print(f"\nCopying head file to: {hds_path}") + shutil.copy2(self.hds_file, hds_path) + + # 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_budget_header_lookup(self): + """ + 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. + """ + result = {} + for rec in self.cbc_obj.recordarray: + result.setdefault((int(rec["kstp"]) - 1, int(rec["kper"]) - 1), rec) + return result + + 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...") + + 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] + totim = float(rec["totim"]) + pertim = float(rec["pertim"]) + delt = float(rec["delt"]) + + if verbose: + print(f" Processing time step {kstpkper}...") + + head = self.hds_obj.get_data(kstpkper=kstpkper) + + 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") + + # 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) + + 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, + } + ) + + # TODO: get_specific_discharge() requires a model object and cannot + # 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. + + sat = get_saturation( + head, self.top, self.botm, self.icelltype_3d, self.hdry, self.hnoflo + ) + + 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 + 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, + } + ) + + 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")" + ) 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