From 25d66970513a961ba00ff29fc549e5854509b33c Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Sat, 12 Nov 2022 08:41:23 -0500 Subject: [PATCH] refactor(utils): move utils from modflow-devtools * add head and budget file writing util fns * add disu creation helper fn * add uniform flow field fn * add head/budget/concentration comparison utility fns * add namefile-reading utility fn * lint tests and notebook utils * add test for HeadUFile.get_ts() * add minimal tests for relocated utils --- autotest/test_binaryfile.py | 44 +- autotest/test_compare.py | 69 ++ autotest/test_gridutil.py | 76 +- autotest/test_mfreadnam.py | 41 + flopy/utils/binaryfile.py | 144 ++++ flopy/utils/compare.py | 1522 +++++++++++++++++++++++++++++++++++ flopy/utils/gridutil.py | 153 ++++ flopy/utils/mfreadnam.py | 490 +++++++++++ 8 files changed, 2537 insertions(+), 2 deletions(-) create mode 100644 autotest/test_compare.py create mode 100644 autotest/test_mfreadnam.py create mode 100644 flopy/utils/compare.py diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index 0cd85c0521..68c305e7ad 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -15,7 +15,12 @@ HeadUFile, Util2d, ) -from flopy.utils.binaryfile import get_headfile_precision +from flopy.utils.binaryfile import ( + get_headfile_precision, + write_budget, + write_head, +) +from flopy.utils.gridutil import uniform_flow_field @pytest.fixture @@ -241,3 +246,40 @@ def test_budgetfile_detect_precision_single(path): def test_budgetfile_detect_precision_double(path): file = CellBudgetFile(path, precision="auto") assert file.realtype == np.float64 + + +def test_write_head(function_tmpdir): + file_path = function_tmpdir / "headfile" + head_data = np.random.random((10, 10)) + + write_head(file_path, head_data) + + assert file_path.is_file() + content = np.fromfile(file_path) + assert np.array_equal(head_data.ravel(), content) + + # TODO: what else needs to be checked here? + + +def test_write_budget(function_tmpdir): + file_path = function_tmpdir / "budgetfile" + + nlay = 3 + nrow = 3 + ncol = 3 + qx = 1.0 + qy = 0.0 + qz = 0.0 + shape = (nlay, nrow, ncol) + spdis, flowja = uniform_flow_field(qx, qy, qz, shape) + + write_budget(file_path, flowja, kstp=0) + assert file_path.is_file() + content1 = np.fromfile(file_path) + + write_budget(file_path, flowja, kstp=1, kper=1, text="text") + assert file_path.is_file() + content2 = np.fromfile(file_path) + + # TODO: why are these the same? + assert np.array_equal(content1, content2) diff --git a/autotest/test_compare.py b/autotest/test_compare.py new file mode 100644 index 0000000000..1e33bc9592 --- /dev/null +++ b/autotest/test_compare.py @@ -0,0 +1,69 @@ +import numpy as np +import pytest + +from flopy.mf6.utils import MfGrdFile +from flopy.utils.compare import _diffmax, _difftol + + +def test_diffmax(): + a1 = np.array([1, 2, 3]) + a2 = np.array([4, 5, 7]) + d, indices = _diffmax(a1, a2) + indices = indices[ + 0 + ] # return value is a tuple of arrays (1 for each dimension) + assert d == 4 + assert list(indices) == [2] + + +def test_difftol(): + a1 = np.array([1, 2, 3]) + a2 = np.array([3, 5, 7]) + d, indices = _difftol(a1, a2, 2.5) + indices = indices[ + 0 + ] # return value is a tuple of arrays (1 for each dimension) + assert d == 4 + print(d, indices) + assert list(indices) == [1, 2] + + +@pytest.mark.skip(reason="todo") +def test_eval_bud_diff(example_data_path): + # get ia from grdfile + mfgrd_test_path = example_data_path / "mfgrd_test" + grb_path = mfgrd_test_path / "nwtp3.dis.grb" + grb = MfGrdFile(str(grb_path), verbose=True) + ia = grb._datadict["IA"] - 1 + + # TODO: create/run minimal model, then compare budget files + + +@pytest.mark.skip(reason="todo") +def test_compare_budget(): + pass + + +@pytest.mark.skip(reason="todo") +def test_compare_swrbudget(): + pass + + +@pytest.mark.skip(reason="todo") +def test_compare_heads(): + pass + + +@pytest.mark.skip(reason="todo") +def test_compare_concs(): + pass + + +@pytest.mark.skip(reason="todo") +def test_compare_stages(): + pass + + +@pytest.mark.skip(reason="todo") +def test_compare(): + pass diff --git a/autotest/test_gridutil.py b/autotest/test_gridutil.py index 7c9640a060..3b66cb4fd7 100644 --- a/autotest/test_gridutil.py +++ b/autotest/test_gridutil.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from flopy.utils.gridutil import get_lni +from flopy.utils.gridutil import get_disu_kwargs, get_lni, uniform_flow_field @pytest.mark.parametrize( @@ -54,3 +54,77 @@ def test_get_lni_infers_layer_count_when_int_ncpl(ncpl, nodes, expected): assert isinstance(lni, list) for i, ln in enumerate(lni): assert ln == expected[i] + + +@pytest.mark.parametrize( + "nlay, nrow, ncol, delr, delc, tp, botm", + [ + ( + 1, + 61, + 61, + np.array(61 * [50]), + np.array(61 * [50]), + np.array([-10]), + np.array([-30.0, -50.0]), + ), + ( + 2, + 61, + 61, + np.array(61 * [50]), + np.array(61 * [50]), + np.array([-10]), + np.array([-30.0, -50.0]), + ), + ], +) +def test_get_disu_kwargs(nlay, nrow, ncol, delr, delc, tp, botm): + kwargs = get_disu_kwargs( + nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, tp=tp, botm=botm + ) + + from pprint import pprint + + pprint(kwargs["area"]) + + assert kwargs["nodes"] == nlay * nrow * ncol + assert kwargs["nvert"] == None + + area = np.array([dr * dc for (dr, dc) in product(delr, delc)], dtype=float) + area = np.array(nlay * [area]).flatten() + assert np.array_equal(kwargs["area"], area) + + # TODO: test other properties + # print(kwargs["iac"]) + # print(kwargs["ihc"]) + # print(kwargs["ja"]) + # print(kwargs["nja"]) + + +@pytest.mark.parametrize( + "qx, qy, qz, nlay, nrow, ncol", + [ + (1, 0, 0, 1, 1, 10), + (0, 1, 0, 1, 1, 10), + (0, 0, 1, 1, 1, 10), + (1, 0, 0, 1, 10, 10), + (1, 0, 0, 2, 10, 10), + (1, 1, 0, 2, 10, 10), + (1, 1, 1, 2, 10, 10), + (2, 1, 1, 2, 10, 10), + ], +) +def test_uniform_flow_field(qx, qy, qz, nlay, nrow, ncol): + shape = nlay, nrow, ncol + spdis, flowja = uniform_flow_field(qx, qy, qz, shape) + + assert spdis.shape == (nlay * nrow * ncol,) + for i, t in enumerate(spdis.flatten()): + assert t[0] == t[1] == i + assert t[3] == qx + assert t[4] == qy + assert t[5] == qz + + # TODO: check flowja + # print(flowja.shape) diff --git a/autotest/test_mfreadnam.py b/autotest/test_mfreadnam.py new file mode 100644 index 0000000000..412b2423a4 --- /dev/null +++ b/autotest/test_mfreadnam.py @@ -0,0 +1,41 @@ +import pytest +from autotest.conftest import get_example_data_path + +from flopy.utils.mfreadnam import get_entries_from_namefile + +_example_data_path = get_example_data_path() + + +@pytest.mark.parametrize( + "path", + [ + _example_data_path / "mf6" / "test001a_Tharmonic" / "mfsim.nam", + _example_data_path / "mf6" / "test001e_UZF_3lay" / "mfsim.nam", + _example_data_path / "mf6-freyberg" / "mfsim.nam", + ], +) +def test_get_entries_from_namefile_mf6(path): + package = "IMS6" + entries = get_entries_from_namefile(path, ftype=package) + assert len(entries) == 1 + + entry = entries[0] + assert path.parent.name in entry[0] + assert entry[1] == package + + +@pytest.mark.skip(reason="only supports mf6 namefiles") +@pytest.mark.parametrize( + "path", + [ + _example_data_path / "mf6-freyberg" / "freyberg.nam", + ], +) +def test_get_entries_from_namefile_mf2005(path): + package = "IC6" + entries = get_entries_from_namefile(path, ftype=package) + assert len(entries) == 1 + + entry = entries[0] + assert path.parent.name in entry[0] + assert entry[1] == package diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index fae2190fa6..aeeeab464d 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile.py @@ -16,6 +16,150 @@ from .gridutil import get_lni +def write_head( + fbin, + data, + kstp=1, + kper=1, + pertim=1.0, + totim=1.0, + text=" HEAD", + ilay=1, +): + dt = np.dtype( + [ + ("kstp", np.int32), + ("kper", np.int32), + ("pertim", np.float64), + ("totim", np.float64), + ("text", "S16"), + ("ncol", np.int32), + ("nrow", np.int32), + ("ilay", np.int32), + ] + ) + nrow = data.shape[0] + ncol = data.shape[1] + h = np.array((kstp, kper, pertim, totim, text, ncol, nrow, ilay), dtype=dt) + h.tofile(fbin) + data.tofile(fbin) + + +def write_budget( + fbin, + data, + kstp=1, + kper=1, + text=" FLOW-JA-FACE", + imeth=1, + delt=1.0, + pertim=1.0, + totim=1.0, + text1id1=" GWF-1", + text2id1=" GWF-1", + text1id2=" GWF-1", + text2id2=" NPF", +): + dt = np.dtype( + [ + ("kstp", np.int32), + ("kper", np.int32), + ("text", "S16"), + ("ndim1", np.int32), + ("ndim2", np.int32), + ("ndim3", np.int32), + ("imeth", np.int32), + ("delt", np.float64), + ("pertim", np.float64), + ("totim", np.float64), + ] + ) + + if imeth == 1: + ndim1 = data.shape[0] + ndim2 = 1 + ndim3 = -1 + h = np.array( + ( + kstp, + kper, + text, + ndim1, + ndim2, + ndim3, + imeth, + delt, + pertim, + totim, + ), + dtype=dt, + ) + h.tofile(fbin) + data.tofile(fbin) + + elif imeth == 6: + ndim1 = 1 + ndim2 = 1 + ndim3 = -1 + h = np.array( + ( + kstp, + kper, + text, + ndim1, + ndim2, + ndim3, + imeth, + delt, + pertim, + totim, + ), + dtype=dt, + ) + h.tofile(fbin) + + # write text1id1, ... + dt = np.dtype( + [ + ("text1id1", "S16"), + ("text1id2", "S16"), + ("text2id1", "S16"), + ("text2id2", "S16"), + ] + ) + h = np.array((text1id1, text1id2, text2id1, text2id2), dtype=dt) + h.tofile(fbin) + + # write ndat (number of floating point columns) + colnames = data.dtype.names + ndat = len(colnames) - 2 + dt = np.dtype([("ndat", np.int32)]) + h = np.array([(ndat,)], dtype=dt) + h.tofile(fbin) + + # write auxiliary column names + naux = ndat - 1 + if naux > 0: + auxtxt = [f"{colname:16}" for colname in colnames[3:]] + auxtxt = tuple(auxtxt) + dt = np.dtype([(colname, "S16") for colname in colnames[3:]]) + h = np.array(auxtxt, dtype=dt) + h.tofile(fbin) + + # write nlist + nlist = data.shape[0] + dt = np.dtype([("nlist", np.int32)]) + h = np.array([(nlist,)], dtype=dt) + h.tofile(fbin) + + # write the data + data.tofile(fbin) + + pass + else: + raise Exception(f"unknown method code {imeth}") + + class BinaryHeader(Header): """ The binary_header class is a class to create headers for MODFLOW diff --git a/flopy/utils/compare.py b/flopy/utils/compare.py new file mode 100644 index 0000000000..c64b98e6d3 --- /dev/null +++ b/flopy/utils/compare.py @@ -0,0 +1,1522 @@ +import os +import textwrap + +import numpy as np + +from flopy.utils.mfreadnam import get_entries_from_namefile + + +def _diffmax(v1, v2): + """Calculate the maximum difference between two vectors. + + Parameters + ---------- + v1 : numpy.ndarray + array of base model results + v2 : numpy.ndarray + array of comparison model results + + Returns + ------- + diffmax : float + absolute value of the maximum difference in v1 and v2 array values + indices : numpy.ndarry + indices where the absolute value of the difference is equal to the + absolute value of the maximum difference. + + """ + if v1.ndim > 1 or v2.ndim > 1: + v1 = v1.flatten() + v2 = v2.flatten() + if v1.size != v2.size: + err = ( + f"Error: calculate_difference v1 size ({v1.size}) " + + f"is not equal to v2 size ({v2.size})" + ) + raise Exception(err) + + diff = abs(v1 - v2) + diffmax = diff.max() + return diffmax, np.where(diff == diffmax) + + +def _difftol(v1, v2, tol): + """Calculate the difference between two arrays relative to a tolerance. + + Parameters + ---------- + v1 : numpy.ndarray + array of base model results + v2 : numpy.ndarray + array of comparison model results + tol : float + tolerance used to evaluate base and comparison models + + Returns + ------- + diffmax : float + absolute value of the maximum difference in v1 and v2 array values + indices : numpy.ndarry + indices where the absolute value of the difference exceed the + specified tolerance. + + """ + if v1.ndim > 1 or v2.ndim > 1: + v1 = v1.flatten() + v2 = v2.flatten() + if v1.size != v2.size: + err = ( + f"Error: calculate_difference v1 size ({v1.size}) " + + f"is not equal to v2 size ({v2.size})" + ) + raise Exception(err) + + diff = abs(v1 - v2) + return diff.max(), np.where(diff > tol) + + +def compare_budget( + namefile1, + namefile2, + max_cumpd=0.01, + max_incpd=0.01, + outfile=None, + files1=None, + files2=None, +): + """Compare the budget results from two simulations. + + Parameters + ---------- + namefile1 : str + namefile path for base model + namefile2 : str + namefile path for comparison model + max_cumpd : float + maximum percent discrepancy allowed for cumulative budget terms + (default is 0.01) + max_incpd : float + maximum percent discrepancy allowed for incremental budget terms + (default is 0.01) + outfile : str + budget comparison output file name. If outfile is None, no + comparison output is saved. (default is None) + files1 : str + base model output file. If files1 is not None, results + will be extracted from files1 and namefile1 will not be used. + (default is None) + files2 : str + comparison model output file. If files2 is not None, results + will be extracted from files2 and namefile2 will not be used. + (default is None) + + Returns + ------- + success : bool + boolean indicating if the difference between budgets are less + than max_cumpd and max_incpd + + """ + try: + import flopy + except: + msg = "flopy not available - cannot use compare_budget" + raise ValueError(msg) + + # headers + headers = ("INCREMENTAL", "CUMULATIVE") + direction = ("IN", "OUT") + + # Get name of list files + lst_file1 = None + if files1 is None: + lst_file = get_entries_from_namefile(namefile1, "list") + lst_file1 = lst_file[0][0] + else: + if isinstance(files1, str): + files1 = [files1] + for file in files1: + if ( + "list" in os.path.basename(file).lower() + or "lst" in os.path.basename(file).lower() + ): + lst_file1 = file + break + lst_file2 = None + if files2 is None: + lst_file = get_entries_from_namefile(namefile2, "list") + lst_file2 = lst_file[0][0] + else: + if isinstance(files2, str): + files2 = [files2] + for file in files2: + if ( + "list" in os.path.basename(file).lower() + or "lst" in os.path.basename(file).lower() + ): + lst_file2 = file + break + # Determine if there are two files to compare + if lst_file1 is None or lst_file2 is None: + print("lst_file1 or lst_file2 is None") + print(f"lst_file1: {lst_file1}") + print(f"lst_file2: {lst_file2}") + return True + + # Open output file + if outfile is not None: + f = open(outfile, "w") + f.write("Created by pymake.autotest.compare\n") + + # Initialize SWR budget objects + lst1obj = flopy.utils.MfusgListBudget(lst_file1) + lst2obj = flopy.utils.MfusgListBudget(lst_file2) + + # Determine if there any SWR entries in the budget file + if not lst1obj.isvalid() or not lst2obj.isvalid(): + return True + + # Get numpy budget tables for lst_file1 + lst1 = [] + lst1.append(lst1obj.get_incremental()) + lst1.append(lst1obj.get_cumulative()) + + # Get numpy budget tables for lst_file2 + lst2 = [] + lst2.append(lst2obj.get_incremental()) + lst2.append(lst2obj.get_cumulative()) + + icnt = 0 + v0 = np.zeros(2, dtype=float) + v1 = np.zeros(2, dtype=float) + err = np.zeros(2, dtype=float) + + # Process cumulative and incremental + for idx in range(2): + if idx > 0: + max_pd = max_cumpd + else: + max_pd = max_incpd + kper = lst1[idx]["stress_period"] + kstp = lst1[idx]["time_step"] + + # Process each time step + for jdx in range(kper.shape[0]): + + err[:] = 0.0 + t0 = lst1[idx][jdx] + t1 = lst2[idx][jdx] + + if outfile is not None: + + maxcolname = 0 + for colname in t0.dtype.names: + maxcolname = max(maxcolname, len(colname)) + + s = 2 * "\n" + s += ( + f"STRESS PERIOD: {kper[jdx] + 1} " + + f"TIME STEP: {kstp[jdx] + 1}" + ) + f.write(s) + + if idx == 0: + f.write("\nINCREMENTAL BUDGET\n") + else: + f.write("\nCUMULATIVE BUDGET\n") + + for i, colname in enumerate(t0.dtype.names): + if i == 0: + s = ( + f"{'Budget Entry':<21} {'Model 1':>21} " + + f"{'Model 2':>21} {'Difference':>21}\n" + ) + f.write(s) + s = 87 * "-" + "\n" + f.write(s) + diff = t0[colname] - t1[colname] + s = ( + f"{colname:<21} {t0[colname]:>21} " + + f"{t1[colname]:>21} {diff:>21}\n" + ) + f.write(s) + + v0[0] = t0["TOTAL_IN"] + v1[0] = t1["TOTAL_IN"] + if v0[0] > 0.0: + err[0] = 100.0 * (v1[0] - v0[0]) / v0[0] + v0[1] = t0["TOTAL_OUT"] + v1[1] = t1["TOTAL_OUT"] + if v0[1] > 0.0: + err[1] = 100.0 * (v1[1] - v0[1]) / v0[1] + for kdx, t in enumerate(err): + if abs(t) > max_pd: + icnt += 1 + if outfile is not None: + e = ( + f'"{headers[idx]} {direction[kdx]}" ' + + f"percent difference ({t})" + + f" for stress period {kper[jdx] + 1} " + + f"and time step {kstp[jdx] + 1} > {max_pd}." + + f" Reference value = {v0[kdx]}. " + + f"Simulated value = {v1[kdx]}." + ) + e = textwrap.fill( + e, + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + f.write(f"{e}\n") + f.write("\n") + + # Close output file + if outfile is not None: + f.close() + + # test for failure + success = True + if icnt > 0: + success = False + return success + + +def compare_swrbudget( + namefile1, + namefile2, + max_cumpd=0.01, + max_incpd=0.01, + outfile=None, + files1=None, + files2=None, +): + """Compare the SWR budget results from two simulations. + + Parameters + ---------- + namefile1 : str + namefile path for base model + namefile2 : str + namefile path for comparison model + max_cumpd : float + maximum percent discrepancy allowed for cumulative budget terms + (default is 0.01) + max_incpd : float + maximum percent discrepancy allowed for incremental budget terms + (default is 0.01) + outfile : str + budget comparison output file name. If outfile is None, no + comparison output is saved. (default is None) + files1 : str + base model output file. If files1 is not None, results + will be extracted from files1 and namefile1 will not be used. + (default is None) + files2 : str + comparison model output file. If files2 is not None, results + will be extracted from files2 and namefile2 will not be used. + (default is None) + + Returns + ------- + success : bool + boolean indicating if the difference between budgets are less + than max_cumpd and max_incpd + + """ + try: + import flopy + except: + msg = "flopy not available - cannot use compare_swrbudget" + raise ValueError(msg) + + # headers + headers = ("INCREMENTAL", "CUMULATIVE") + direction = ("IN", "OUT") + + # Get name of list files + list1 = None + if files1 is None: + lst = get_entries_from_namefile(namefile1, "list") + list1 = lst[0][0] + else: + for file in files1: + if ( + "list" in os.path.basename(file).lower() + or "lst" in os.path.basename(file).lower() + ): + list1 = file + break + list2 = None + if files2 is None: + lst = get_entries_from_namefile(namefile2, "list") + list2 = lst[0][0] + else: + for file in files2: + if ( + "list" in os.path.basename(file).lower() + or "lst" in os.path.basename(file).lower() + ): + list2 = file + break + # Determine if there are two files to compare + if list1 is None or list2 is None: + return True + + # Initialize SWR budget objects + lst1obj = flopy.utils.SwrListBudget(list1) + lst2obj = flopy.utils.SwrListBudget(list2) + + # Determine if there any SWR entries in the budget file + if not lst1obj.isvalid() or not lst2obj.isvalid(): + return True + + # Get numpy budget tables for list1 + lst1 = [] + lst1.append(lst1obj.get_incremental()) + lst1.append(lst1obj.get_cumulative()) + + # Get numpy budget tables for list2 + lst2 = [] + lst2.append(lst2obj.get_incremental()) + lst2.append(lst2obj.get_cumulative()) + + icnt = 0 + v0 = np.zeros(2, dtype=float) + v1 = np.zeros(2, dtype=float) + err = np.zeros(2, dtype=float) + + # Open output file + if outfile is not None: + f = open(outfile, "w") + f.write("Created by pymake.autotest.compare\n") + + # Process cumulative and incremental + for idx in range(2): + if idx > 0: + max_pd = max_cumpd + else: + max_pd = max_incpd + kper = lst1[idx]["stress_period"] + kstp = lst1[idx]["time_step"] + + # Process each time step + for jdx in range(kper.shape[0]): + + err[:] = 0.0 + t0 = lst1[idx][jdx] + t1 = lst2[idx][jdx] + + if outfile is not None: + + maxcolname = 0 + for colname in t0.dtype.names: + maxcolname = max(maxcolname, len(colname)) + + s = 2 * "\n" + s += ( + f"STRESS PERIOD: {kper[jdx] + 1} " + + f"TIME STEP: {kstp[jdx] + 1}" + ) + f.write(s) + + if idx == 0: + f.write("\nINCREMENTAL BUDGET\n") + else: + f.write("\nCUMULATIVE BUDGET\n") + + for i, colname in enumerate(t0.dtype.names): + if i == 0: + s = ( + f"{'Budget Entry':<21} {'Model 1':>21} " + + f"{'Model 2':>21} {'Difference':>21}\n" + ) + f.write(s) + s = 87 * "-" + "\n" + f.write(s) + diff = t0[colname] - t1[colname] + s = ( + f"{colname:<21} {t0[colname]:>21} " + + f"{t1[colname]:>21} {diff:>21}\n" + ) + f.write(s) + + v0[0] = t0["TOTAL_IN"] + v1[0] = t1["TOTAL_IN"] + if v0[0] > 0.0: + err[0] = 100.0 * (v1[0] - v0[0]) / v0[0] + v0[1] = t0["TOTAL_OUT"] + v1[1] = t1["TOTAL_OUT"] + if v0[1] > 0.0: + err[1] = 100.0 * (v1[1] - v0[1]) / v0[1] + for kdx, t in enumerate(err): + if abs(t) > max_pd: + icnt += 1 + e = ( + f'"{headers[idx]} {direction[kdx]}" ' + + f"percent difference ({t})" + + f" for stress period {kper[jdx] + 1} " + + f"and time step {kstp[jdx] + 1} > {max_pd}." + + f" Reference value = {v0[kdx]}. " + + f"Simulated value = {v1[kdx]}." + ) + e = textwrap.fill( + e, + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + f.write(f"{e}\n") + f.write("\n") + + # Close output file + if outfile is not None: + f.close() + + # test for failure + success = True + if icnt > 0: + success = False + return success + + +def compare_heads( + namefile1, + namefile2, + precision="auto", + text="head", + text2=None, + htol=0.001, + outfile=None, + files1=None, + files2=None, + difftol=False, + verbose=False, + exfile=None, + exarr=None, + maxerr=None, +): + """Compare the head results from two simulations. + + Parameters + ---------- + namefile1 : str + namefile path for base model + namefile2 : str + namefile path for comparison model + precision : str + precision for binary head file ("auto", "single", or "double") + default is "auto" + htol : float + maximum allowed head difference (default is 0.001) + outfile : str + head comparison output file name. If outfile is None, no + comparison output is saved. (default is None) + files1 : str + base model output file. If files1 is not None, results + will be extracted from files1 and namefile1 will not be used. + (default is None) + files2 : str + comparison model output file. If files2 is not None, results + will be extracted from files2 and namefile2 will not be used. + (default is None) + difftol : bool + boolean determining if the absolute value of the head + difference greater than htol should be evaluated (default is False) + verbose : bool + boolean indicating if verbose output should be written to the + terminal (default is False) + exfile : str + path to a file with exclusion array data. Head differences will not + be evaluated where exclusion array values are greater than zero. + (default is None) + exarr : numpy.ndarry + exclusion array. Head differences will not be evaluated where + exclusion array values are greater than zero. (default is None). + maxerr : int + maximum number of head difference greater than htol that should be + reported. If maxerr is None, all head difference greater than htol + will be reported. (default is None) + + Returns + ------- + success : bool + boolean indicating if the head differences are less than htol. + + """ + try: + import flopy + except: + msg = "flopy not available - cannot use compare_heads" + raise ValueError(msg) + + if text2 is None: + text2 = text + + dbs = "DATA(BINARY)" + + # Get head info for namefile1 + hfpth1 = None + status1 = dbs + if files1 is None: + # Get oc info, and return if OC not included in models + ocf1 = get_entries_from_namefile(namefile1, "OC") + if ocf1[0][0] is None: + return True + + hu1, hfpth1, du1, _ = flopy.modflow.ModflowOc.get_ocoutput_units( + ocf1[0][0] + ) + if text.lower() == "head": + iut = hu1 + elif text.lower() == "drawdown": + iut = du1 + if iut != 0: + entries = get_entries_from_namefile(namefile1, unit=abs(iut)) + hfpth1, status1 = entries[0][0], entries[0][1] + + else: + if isinstance(files1, str): + files1 = [files1] + for file in files1: + if text.lower() == "head": + if ( + "hds" in os.path.basename(file).lower() + or "hed" in os.path.basename(file).lower() + ): + hfpth1 = file + break + elif text.lower() == "drawdown": + if "ddn" in os.path.basename(file).lower(): + hfpth1 = file + break + elif text.lower() == "concentration": + if "ucn" in os.path.basename(file).lower(): + hfpth1 = file + break + else: + hfpth1 = file + break + + # Get head info for namefile2 + hfpth2 = None + status2 = dbs + if files2 is None: + # Get oc info, and return if OC not included in models + ocf2 = get_entries_from_namefile(namefile2, "OC") + if ocf2[0][0] is None: + return True + + hu2, hfpth2, du2, dfpth2 = flopy.modflow.ModflowOc.get_ocoutput_units( + ocf2[0][0] + ) + if text.lower() == "head": + iut = hu2 + elif text.lower() == "drawdown": + iut = du2 + if iut != 0: + entries = get_entries_from_namefile(namefile2, unit=abs(iut)) + hfpth2, status2 = entries[0][0], entries[0][1] + else: + if isinstance(files2, str): + files2 = [files2] + for file in files2: + if text2.lower() == "head": + if ( + "hds" in os.path.basename(file).lower() + or "hed" in os.path.basename(file).lower() + ): + hfpth2 = file + break + elif text2.lower() == "drawdown": + if "ddn" in os.path.basename(file).lower(): + hfpth2 = file + break + elif text2.lower() == "concentration": + if "ucn" in os.path.basename(file).lower(): + hfpth2 = file + break + else: + hfpth2 = file + break + + # confirm that there are two files to compare + if hfpth1 is None or hfpth2 is None: + print("hfpth1 or hfpth2 is None") + print(f"hfpth1: {hfpth1}") + print(f"hfpth2: {hfpth2}") + return True + + # make sure the file paths exist + if not os.path.isfile(hfpth1) or not os.path.isfile(hfpth2): + print("hfpth1 or hfpth2 is not a file") + print(f"hfpth1 isfile: {os.path.isfile(hfpth1)}") + print(f"hfpth2 isfile: {os.path.isfile(hfpth2)}") + return False + + # Open output file + if outfile is not None: + f = open(outfile, "w") + f.write("Created by pymake.autotest.compare\n") + f.write(f"Performing {text.upper()} to {text2.upper()} comparison\n") + + if exfile is not None: + f.write(f"Using exclusion file {exfile}\n") + if exarr is not None: + f.write("Using exclusion array\n") + + msg = f"{hfpth1} is a " + if status1 == dbs: + msg += "binary file." + else: + msg += "ascii file." + f.write(msg + "\n") + msg = f"{hfpth2} is a " + if status2 == dbs: + msg += "binary file." + else: + msg += "ascii file." + f.write(msg + "\n") + + # Process exclusion data + exd = None + # get data from exclusion file + if exfile is not None: + e = None + if isinstance(exfile, str): + try: + exd = np.genfromtxt(exfile).flatten() + except: + e = ( + "Could not read exclusion " + + f"file {os.path.basename(exfile)}" + ) + print(e) + return False + else: + e = "exfile is not a valid file path" + print(e) + return False + + # process exclusion array + if exarr is not None: + e = None + if isinstance(exarr, np.ndarray): + if exd is None: + exd = exarr.flatten() + else: + exd += exarr.flatten() + else: + e = "exarr is not a numpy array" + print(e) + return False + + # Get head objects + status1 = status1.upper() + unstructured1 = False + if status1 == dbs: + headobj1 = flopy.utils.HeadFile( + hfpth1, precision=precision, verbose=verbose, text=text + ) + txt = headobj1.recordarray["text"][0] + if isinstance(txt, bytes): + txt = txt.decode("utf-8") + if "HEADU" in txt: + unstructured1 = True + headobj1 = flopy.utils.HeadUFile( + hfpth1, precision=precision, verbose=verbose + ) + else: + headobj1 = flopy.utils.FormattedHeadFile( + hfpth1, verbose=verbose, text=text + ) + + status2 = status2.upper() + unstructured2 = False + if status2 == dbs: + headobj2 = flopy.utils.HeadFile( + hfpth2, precision=precision, verbose=verbose, text=text2 + ) + txt = headobj2.recordarray["text"][0] + if isinstance(txt, bytes): + txt = txt.decode("utf-8") + if "HEADU" in txt: + unstructured2 = True + headobj2 = flopy.utils.HeadUFile( + hfpth2, precision=precision, verbose=verbose + ) + else: + headobj2 = flopy.utils.FormattedHeadFile( + hfpth2, verbose=verbose, text=text2 + ) + + # get times + times1 = headobj1.get_times() + times2 = headobj2.get_times() + for (t1, t2) in zip(times1, times2): + if not np.allclose([t1], [t2]): + msg = "times in two head files are not " + f"equal ({t1},{t2})" + raise ValueError(msg) + + kstpkper = headobj1.get_kstpkper() + + line_separator = 15 * "-" + header = ( + f"{' ':>15s} {' ':>15s} {'MAXIMUM':>15s} {'EXCEEDS':>15s}\n" + + f"{'STRESS PERIOD':>15s} {'TIME STEP':>15s} " + + f"{'HEAD DIFFERENCE':>15s} {'CRITERIA':>15s}\n" + + f"{line_separator:>15s} {line_separator:>15s} " + + f"{line_separator:>15s} {line_separator:>15s}\n" + ) + + if verbose: + print(f"Comparing results for {len(times1)} times") + + icnt = 0 + # Process cumulative and incremental + for idx, (t1, t2) in enumerate(zip(times1, times2)): + h1 = headobj1.get_data(totim=t1) + if unstructured1: + temp = np.array([]) + for a in h1: + temp = np.hstack((temp, a)) + h1 = temp + h2 = headobj2.get_data(totim=t2) + if unstructured2: + temp = np.array([]) + for a in h2: + temp = np.hstack((temp, a)) + h2 = temp + + if exd is not None: + # reshape exd to the shape of the head arrays + if idx == 0: + e = ( + f"shape of exclusion data ({exd.shape})" + + "can not be reshaped to the size of the " + + f"head arrays ({h1.shape})" + ) + if h1.flatten().shape != exd.shape: + raise ValueError(e) + exd = exd.reshape(h1.shape) + iexd = exd > 0 + + # reset h1 and h2 to the same value in the excluded area + h1[iexd] = 0.0 + h2[iexd] = 0.0 + + if difftol: + diffmax, indices = _difftol(h1, h2, htol) + else: + diffmax, indices = _diffmax(h1, h2) + + if outfile is not None: + if idx < 1: + f.write(header) + if diffmax > htol: + sexceed = "*" + else: + sexceed = "" + kk1 = kstpkper[idx][1] + 1 + kk0 = kstpkper[idx][0] + 1 + f.write(f"{kk1:15d} {kk0:15d} {diffmax:15.6g} {sexceed:15s}\n") + + if diffmax >= htol: + icnt += 1 + if outfile is not None: + if difftol: + ee = ( + "Maximum absolute head difference " + + f"({diffmax}) -- " + + f"{htol} tolerance exceeded at " + + f"{indices[0].shape[0]} node location(s)" + ) + else: + ee = ( + "Maximum absolute head difference " + + f"({diffmax}) exceeded " + + f"at {indices[0].shape[0]} node location(s)" + ) + e = textwrap.fill( + ee + ":", + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + + if verbose: + f.write(f"{ee}\n") + print(ee + f" at time {t1}") + + e = "" + ncells = h1.flatten().shape[0] + fmtn = "{:" + f"{len(str(ncells))}" + "d}" + for itupe in indices: + for jdx, ind in enumerate(itupe): + iv = np.unravel_index(ind, h1.shape) + iv = tuple(i + 1 for i in iv) + v1 = h1.flatten()[ind] + v2 = h2.flatten()[ind] + d12 = v1 - v2 + # e += ' ' + fmtn.format(jdx + 1) + ' node: ' + # e += fmtn.format(ind + 1) # convert to one-based + e += " " + fmtn.format(jdx + 1) + e += f" {iv}" + e += " -- " + e += f"h1: {v1:20} " + e += f"h2: {v2:20} " + e += f"diff: {d12:20}\n" + if isinstance(maxerr, int): + if jdx + 1 >= maxerr: + break + if verbose: + f.write(f"{e}\n") + # Write header again, unless it is the last record + if verbose: + if idx + 1 < len(times1): + f.write(f"\n{header}") + + # Close output file + if outfile is not None: + f.close() + + # test for failure + success = True + if icnt > 0: + success = False + return success + + +def compare_concentrations( + namefile1, + namefile2, + precision="auto", + ctol=0.001, + outfile=None, + files1=None, + files2=None, + difftol=False, + verbose=False, +): + """Compare the mt3dms and mt3dusgs concentration results from two + simulations. + + Parameters + ---------- + namefile1 : str + namefile path for base model + namefile2 : str + namefile path for comparison model + precision : str + precision for binary head file ("auto", "single", or "double") + default is "auto" + ctol : float + maximum allowed concentration difference (default is 0.001) + outfile : str + concentration comparison output file name. If outfile is None, no + comparison output is saved. (default is None) + files1 : str + base model output file. If files1 is not None, results + will be extracted from files1 and namefile1 will not be used. + (default is None) + files2 : str + comparison model output file. If files2 is not None, results + will be extracted from files2 and namefile2 will not be used. + (default is None) + difftol : bool + boolean determining if the absolute value of the concentration + difference greater than ctol should be evaluated (default is False) + verbose : bool + boolean indicating if verbose output should be written to the + terminal (default is False) + + Returns + ------- + success : bool + boolean indicating if the concentration differences are less than + ctol. + + Returns + ------- + + """ + try: + import flopy + except: + msg = "flopy not available - cannot use compare_concs" + raise ValueError(msg) + + # list of valid extensions + valid_ext = ["ucn"] + + # Get info for first ucn file + ufpth1 = None + if files1 is None: + for ext in valid_ext: + ucn = get_entries_from_namefile(namefile1, extension=ext) + ufpth = ucn[0][0] + if ufpth is not None: + ufpth1 = ufpth + break + if ufpth1 is None: + ufpth1 = os.path.join(os.path.dirname(namefile1), "MT3D001.UCN") + else: + if isinstance(files1, str): + files1 = [files1] + for file in files1: + for ext in valid_ext: + if ext in os.path.basename(file).lower(): + ufpth1 = file + break + + # Get info for second ucn file + ufpth2 = None + if files2 is None: + for ext in valid_ext: + ucn = get_entries_from_namefile(namefile2, extension=ext) + ufpth = ucn[0][0] + if ufpth is not None: + ufpth2 = ufpth + break + if ufpth2 is None: + ufpth2 = os.path.join(os.path.dirname(namefile2), "MT3D001.UCN") + else: + if isinstance(files2, str): + files2 = [files2] + for file in files2: + for ext in valid_ext: + if ext in os.path.basename(file).lower(): + ufpth2 = file + break + + # confirm that there are two files to compare + if ufpth1 is None or ufpth2 is None: + if ufpth1 is None: + print(" UCN file 1 not set") + if ufpth2 is None: + print(" UCN file 2 not set") + return True + + if not os.path.isfile(ufpth1) or not os.path.isfile(ufpth2): + if not os.path.isfile(ufpth1): + print(f" {ufpth1} does not exist") + if not os.path.isfile(ufpth2): + print(f" {ufpth2} does not exist") + return True + + # Open output file + if outfile is not None: + f = open(outfile, "w") + f.write("Created by pymake.autotest.compare_concs\n") + + # Get stage objects + uobj1 = flopy.utils.UcnFile(ufpth1, precision=precision, verbose=verbose) + uobj2 = flopy.utils.UcnFile(ufpth2, precision=precision, verbose=verbose) + + # get times + times1 = uobj1.get_times() + times2 = uobj2.get_times() + nt1 = len(times1) + nt2 = len(times2) + nt = min(nt1, nt2) + + for (t1, t2) in zip(times1, times2): + if not np.allclose([t1], [t2]): + msg = f"times in two ucn files are not equal ({t1},{t2})" + raise ValueError(msg) + + if nt == nt1: + kstpkper = uobj1.get_kstpkper() + else: + kstpkper = uobj2.get_kstpkper() + + line_separator = 15 * "-" + header = ( + f"{' ':>15s} {' ':>15s} {'MAXIMUM':>15s}\n" + + f"{'STRESS PERIOD':>15s} {'TIME STEP':>15s} " + + f"{'CONC DIFFERENCE':>15s}\n" + + f"{line_separator:>15s} " + + f"{line_separator:>15s} " + + f"{line_separator:>15s}\n" + ) + + if verbose: + print(f"Comparing results for {len(times1)} times") + + icnt = 0 + # Process cumulative and incremental + for idx, time in enumerate(times1[0:nt]): + try: + u1 = uobj1.get_data(totim=time) + u2 = uobj2.get_data(totim=time) + + if difftol: + diffmax, indices = _difftol(u1, u2, ctol) + else: + diffmax, indices = _diffmax(u1, u2) + + if outfile is not None: + if idx < 1: + f.write(header) + f.write( + f"{kstpkper[idx][1] + 1:15d} " + + f"{kstpkper[idx][0] + 1:15d} " + + f"{diffmax:15.6g}\n" + ) + + if diffmax >= ctol: + icnt += 1 + if outfile is not None: + if difftol: + ee = ( + f"Maximum concentration difference ({diffmax})" + + f" -- {ctol} tolerance exceeded at " + + f"{indices[0].shape[0]} node location(s)" + ) + else: + ee = ( + "Maximum concentration difference " + + f"({diffmax}) exceeded " + + f"at {indices[0].shape[0]} node location(s)" + ) + e = textwrap.fill( + ee + ":", + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + f.write(f"{e}\n") + if verbose: + print(ee + f" at time {time}") + e = "" + for itupe in indices: + for ind in itupe: + e += f"{ind + 1} " # convert to one-based + e = textwrap.fill( + e, + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + f.write(f"{e}\n") + # Write header again, unless it is the last record + if idx + 1 < len(times1): + f.write(f"\n{header}") + except: + print(f" could not process time={time}") + print(" terminating ucn processing...") + break + + # Close output file + if outfile is not None: + f.close() + + # test for failure + success = True + if icnt > 0: + success = False + return success + + +def compare_stages( + namefile1=None, + namefile2=None, + files1=None, + files2=None, + htol=0.001, + outfile=None, + difftol=False, + verbose=False, +): + """Compare SWR process stage results from two simulations. + + Parameters + ---------- + namefile1 : str + namefile path for base model + namefile2 : str + namefile path for comparison model + precision : str + precision for binary head file ("auto", "single", or "double") + default is "auto" + htol : float + maximum allowed stage difference (default is 0.001) + outfile : str + head comparison output file name. If outfile is None, no + comparison output is saved. (default is None) + files1 : str + base model output file. If files1 is not None, results + will be extracted from files1 and namefile1 will not be used. + (default is None) + files2 : str + comparison model output file. If files2 is not None, results + will be extracted from files2 and namefile2 will not be used. + (default is None) + difftol : bool + boolean determining if the absolute value of the stage + difference greater than htol should be evaluated (default is False) + verbose : bool + boolean indicating if verbose output should be written to the + terminal (default is False) + + Returns + ------- + success : bool + boolean indicating if the stage differences are less than htol. + + """ + try: + import flopy + except: + msg = "flopy not available - cannot use compare_stages" + raise ValueError(msg) + + # list of valid extensions + valid_ext = ["stg"] + + # Get info for first stage file + sfpth1 = None + if namefile1 is not None: + for ext in valid_ext: + stg = get_entries_from_namefile(namefile1, extension=ext) + sfpth = stg[0][0] + if sfpth is not None: + sfpth1 = sfpth + break + elif files1 is not None: + if isinstance(files1, str): + files1 = [files1] + for file in files1: + for ext in valid_ext: + if ext in os.path.basename(file).lower(): + sfpth1 = file + break + + # Get info for second stage file + sfpth2 = None + if namefile2 is not None: + for ext in valid_ext: + stg = get_entries_from_namefile(namefile2, extension=ext) + sfpth = stg[0][0] + if sfpth is not None: + sfpth2 = sfpth + break + elif files2 is not None: + if isinstance(files2, str): + files2 = [files2] + for file in files2: + for ext in valid_ext: + if ext in os.path.basename(file).lower(): + sfpth2 = file + break + + # confirm that there are two files to compare + if sfpth1 is None or sfpth2 is None: + print("spth1 or spth2 is None") + print(f"spth1: {sfpth1}") + print(f"spth2: {sfpth2}") + return False + + if not os.path.isfile(sfpth1) or not os.path.isfile(sfpth2): + print("spth1 or spth2 is not a file") + print(f"spth1 isfile: {os.path.isfile(sfpth1)}") + print(f"spth2 isfile: {os.path.isfile(sfpth2)}") + return False + + # Open output file + if outfile is not None: + f = open(outfile, "w") + f.write("Created by pymake.autotest.compare_stages\n") + + # Get stage objects + sobj1 = flopy.utils.SwrStage(sfpth1, verbose=verbose) + sobj2 = flopy.utils.SwrStage(sfpth2, verbose=verbose) + + # get totim + times1 = sobj1.get_times() + + # get kswr, kstp, and kper + kk = sobj1.get_kswrkstpkper() + + line_separator = 15 * "-" + header = ( + f"{' ':>15s} {' ':>15s} {' ':>15s} {'MAXIMUM':>15s}\n" + + f"{'STRESS PERIOD':>15s} " + + f"{'TIME STEP':>15s} " + + f"{'SWR TIME STEP':>15s} " + + f"{'STAGE DIFFERENCE':>15s}\n" + + f"{line_separator:>15s} " + + f"{line_separator:>15s} " + + f"{line_separator:>15s} " + + f"{line_separator:>15s}\n" + ) + + if verbose: + print(f"Comparing results for {len(times1)} times") + + icnt = 0 + # Process stage data + for idx, (kon, time) in enumerate(zip(kk, times1)): + s1 = sobj1.get_data(totim=time) + s2 = sobj2.get_data(totim=time) + + if s1 is None or s2 is None: + continue + + s1 = s1["stage"] + s2 = s2["stage"] + + if difftol: + diffmax, indices = _difftol(s1, s2, htol) + else: + diffmax, indices = _diffmax(s1, s2) + + if outfile is not None: + if idx < 1: + f.write(header) + f.write( + f"{kon[2] + 1:15d} " + + f"{kon[1] + 1:15d} " + + f"{kon[0] + 1:15d} " + + f"{diffmax:15.6g}\n" + ) + + if diffmax >= htol: + icnt += 1 + if outfile is not None: + if difftol: + ee = ( + f"Maximum head difference ({diffmax}) -- " + + f"{htol} tolerance exceeded at " + + f"{indices[0].shape[0]} node location(s)" + ) + else: + ee = ( + "Maximum head difference " + + f"({diffmax}) exceeded " + + f"at {indices[0].shape[0]} node location(s):" + ) + e = textwrap.fill( + ee + ":", + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + f.write(f"{e}\n") + if verbose: + print(ee + f" at time {time}") + e = "" + for itupe in indices: + for ind in itupe: + e += f"{ind + 1} " # convert to one-based + e = textwrap.fill( + e, + width=70, + initial_indent=" ", + subsequent_indent=" ", + ) + f.write(f"{e}\n") + # Write header again, unless it is the last record + if idx + 1 < len(times1): + f.write(f"\n{header}") + + # Close output file + if outfile is not None: + f.close() + + # test for failure + success = True + if icnt > 0: + success = False + return success + + +def compare( + namefile1, + namefile2, + precision="auto", + max_cumpd=0.01, + max_incpd=0.01, + htol=0.001, + outfile1=None, + outfile2=None, + files1=None, + files2=None, +): + """Compare the budget and head results for two MODFLOW-based model + simulations. + + Parameters + ---------- + namefile1 : str + namefile path for base model + namefile2 : str + namefile path for comparison model + precision : str + precision for binary head file ("auto", "single", or "double") + default is "auto" + max_cumpd : float + maximum percent discrepancy allowed for cumulative budget terms + (default is 0.01) + max_incpd : float + maximum percent discrepancy allowed for incremental budget terms + (default is 0.01) + htol : float + maximum allowed head difference (default is 0.001) + outfile1 : str + budget comparison output file name. If outfile1 is None, no budget + comparison output is saved. (default is None) + outfile2 : str + head comparison output file name. If outfile2 is None, no head + comparison output is saved. (default is None) + files1 : str + base model output file. If files1 is not None, results + will be extracted from files1 and namefile1 will not be used. + (default is None) + files2 : str + comparison model output file. If files2 is not None, results + will be extracted from files2 and namefile2 will not be used. + (default is None) + + Returns + ------- + success : bool + boolean indicating if the budget and head differences are less than + max_cumpd, max_incpd, and htol. + + """ + + # Compare budgets from the list files in namefile1 and namefile2 + success1 = compare_budget( + namefile1, + namefile2, + max_cumpd=max_cumpd, + max_incpd=max_incpd, + outfile=outfile1, + files1=files1, + files2=files2, + ) + success2 = compare_heads( + namefile1, + namefile2, + precision=precision, + htol=htol, + outfile=outfile2, + files1=files1, + files2=files2, + ) + success = False + if success1 and success2: + success = True + return success + + +def eval_bud_diff(fpth, b0, b1, ia=None, dtol=1e-6): + # To use this eval_bud_diff function on a gwf or gwt budget file, + # the function may need ia, in order to exclude comparison of the residual + # term, which is stored in the diagonal position of the flowja array. + # The following code can be used to extract ia from the grb file. + # get ia/ja from binary grid file + # fname = '{}.dis.grb'.format(os.path.basename(sim.name)) + # fpth = os.path.join(sim.simpath, fname) + # grbobj = flopy.mf6.utils.MfGrdFile(fpth) + # ia = grbobj._datadict['IA'] - 1 + + diffmax = 0.0 + difftag = "None" + difftime = None + fail = False + + # build list of cbc data to retrieve + avail = b0.get_unique_record_names() + + # initialize list for storing totals for each budget term terms + cbc_keys = [] + for t in avail: + if isinstance(t, bytes): + t = t.decode() + t = t.strip() + cbc_keys.append(t) + + # open a summary file and write header + f = open(fpth, "w") + line = f"{'Time':15s}" + line += f" {'Datatype':15s}" + line += f" {'File 1':15s}" + line += f" {'File 2':15s}" + line += f" {'Difference':15s}" + f.write(line + "\n") + f.write(len(line) * "-" + "\n") + + # get data from cbc file + kk = b0.get_kstpkper() + times = b0.get_times() + for idx, (k, t) in enumerate(zip(kk, times)): + v0sum = 0.0 + v1sum = 0.0 + for key in cbc_keys: + v0 = b0.get_data(kstpkper=k, text=key)[0] + v1 = b1.get_data(kstpkper=k, text=key)[0] + if isinstance(v0, np.recarray): + v0 = v0["q"].sum() + v1 = v1["q"].sum() + else: + v0 = v0.flatten() + v1 = v1.flatten() + if key == "FLOW-JA-FACE": + # Set residual (stored in diagonal of flowja) to zero + if ia is None: + raise Exception("ia is required for model flowja") + idiagidx = ia[:-1] + v0[idiagidx] = 0.0 + v1[idiagidx] = 0.0 + v0 = v0.sum() + v1 = v1.sum() + + # sum all of the values + if key != "AUXILIARY": + v0sum += v0 + v1sum += v1 + + diff = v0 - v1 + if abs(diff) > abs(diffmax): + diffmax = diff + difftag = key + difftime = t + if abs(diff) > dtol: + fail = True + line = f"{t:15g}" + line += f" {key:15s}" + line += f" {v0:15g}" + line += f" {v1:15g}" + line += f" {diff:15g}" + f.write(line + "\n") + + # evaluate the sums + diff = v0sum - v1sum + if abs(diff) > dtol: + fail = True + line = f"{t:15g}" + line += f" {'TOTAL':15s}" + line += f" {v0sum:15g}" + line += f" {v1sum:15g}" + line += f" {diff:15g}" + f.write(line + "\n") + + msg = f"\nSummary of changes in {os.path.basename(fpth)}\n" + msg += "-" * 72 + "\n" + msg += f"Maximum cbc difference: {diffmax}\n" + msg += f"Maximum cbc difference time: {difftime}\n" + msg += f"Maximum cbc datatype: {difftag}\n" + if fail: + msg += f"Maximum cbc criteria exceeded: {dtol}" + assert not fail, msg + + # close summary file and print the final message + f.close() + print(msg) + + msg = f"sum of first cbc file flows ({v0sum}) " + f"exceeds dtol ({dtol})" + assert abs(v0sum) < dtol, msg + + msg = f"sum of second cbc file flows ({v1sum}) " + f"exceeds dtol ({dtol})" + assert abs(v1sum) < dtol, msg diff --git a/flopy/utils/gridutil.py b/flopy/utils/gridutil.py index 4bb5a109d6..953ee97a68 100644 --- a/flopy/utils/gridutil.py +++ b/flopy/utils/gridutil.py @@ -56,3 +56,156 @@ def get_lni(ncpl, nodes) -> List[Tuple[int, int]]: tuples.append((layer + 1, 0) if correct else (layer, nidx)) return tuples + + +def get_disu_kwargs(nlay, nrow, ncol, delr, delc, tp, botm): + """ + Simple utility for creating args needed to construct + a disu package + """ + + def get_nn(k, i, j): + return k * nrow * ncol + i * ncol + j + + nodes = nlay * nrow * ncol + iac = np.zeros((nodes), dtype=int) + ja = [] + area = np.zeros((nodes), dtype=float) + top = np.zeros((nodes), dtype=float) + bot = np.zeros((nodes), dtype=float) + ihc = [] + cl12 = [] + hwva = [] + for k in range(nlay): + for i in range(nrow): + for j in range(ncol): + # diagonal + n = get_nn(k, i, j) + ja.append(n) + iac[n] += 1 + area[n] = delr[i] * delc[j] + ihc.append(n + 1) + cl12.append(n + 1) + hwva.append(n + 1) + if k == 0: + top[n] = tp + else: + top[n] = botm[k - 1] + bot[n] = botm[k] + # up + if k > 0: + ja.append(get_nn(k - 1, i, j)) + iac[n] += 1 + ihc.append(0) + dz = botm[k - 1] - botm[k] + cl12.append(0.5 * dz) + hwva.append(delr[i] * delc[j]) + # back + if i > 0: + ja.append(get_nn(k, i - 1, j)) + iac[n] += 1 + ihc.append(1) + cl12.append(0.5 * delc[i]) + hwva.append(delr[j]) + # left + if j > 0: + ja.append(get_nn(k, i, j - 1)) + iac[n] += 1 + ihc.append(1) + cl12.append(0.5 * delr[j]) + hwva.append(delc[i]) + # right + if j < ncol - 1: + ja.append(get_nn(k, i, j + 1)) + iac[n] += 1 + ihc.append(1) + cl12.append(0.5 * delr[j]) + hwva.append(delc[i]) + # front + if i < nrow - 1: + ja.append(get_nn(k, i + 1, j)) + iac[n] += 1 + ihc.append(1) + cl12.append(0.5 * delc[i]) + hwva.append(delr[j]) + # bottom + if k < nlay - 1: + ja.append(get_nn(k + 1, i, j)) + iac[n] += 1 + ihc.append(0) + if k == 0: + dz = tp - botm[k] + else: + dz = botm[k - 1] - botm[k] + cl12.append(0.5 * dz) + hwva.append(delr[i] * delc[j]) + ja = np.array(ja, dtype=int) + nja = ja.shape[0] + hwva = np.array(hwva, dtype=float) + kw = {} + kw["nodes"] = nodes + kw["nja"] = nja + kw["nvert"] = None + kw["top"] = top + kw["bot"] = bot + kw["area"] = area + kw["iac"] = iac + kw["ja"] = ja + kw["ihc"] = ihc + kw["cl12"] = cl12 + kw["hwva"] = hwva + return kw + + +def uniform_flow_field(qx, qy, qz, shape, delr=None, delc=None, delv=None): + nlay, nrow, ncol = shape + + # create spdis array for the uniform flow field + dt = np.dtype( + [ + ("ID1", np.int32), + ("ID2", np.int32), + ("FLOW", np.float64), + ("QX", np.float64), + ("QY", np.float64), + ("QZ", np.float64), + ] + ) + spdis = np.array( + [(id1, id1, 0.0, qx, qy, qz) for id1 in range(nlay * nrow * ncol)], + dtype=dt, + ) + + # create the flowja array for the uniform flow field (assume top-bot = 1) + flowja = [] + if delr is None: + delr = 1.0 + if delc is None: + delc = 1.0 + if delv is None: + delv = 1.0 + for k in range(nlay): + for i in range(nrow): + for j in range(ncol): + # diagonal + flowja.append(0.0) + # up + if k > 0: + flowja.append(-qz * delr * delc) + # back + if i > 0: + flowja.append(-qy * delr * delv) + # left + if j > 0: + flowja.append(qx * delc * delv) + # right + if j < ncol - 1: + flowja.append(-qx * delc * delv) + # front + if i < nrow - 1: + flowja.append(qy * delr * delv) + # bottom + if k < nlay - 1: + flowja.append(qz * delr * delc) + flowja = np.array(flowja, dtype=np.float64) + return spdis, flowja diff --git a/flopy/utils/mfreadnam.py b/flopy/utils/mfreadnam.py index 9ac6c07e4a..c50a9497dc 100644 --- a/flopy/utils/mfreadnam.py +++ b/flopy/utils/mfreadnam.py @@ -7,7 +7,10 @@ `_. """ +import os +from os import PathLike from pathlib import Path, PurePosixPath, PureWindowsPath +from typing import List, Tuple class NamData: @@ -266,3 +269,490 @@ def attribs_from_namfile_header(namefile): except: print(f" could not parse start in {namefile}") return defaults + + +def get_entries_from_namefile( + path: PathLike, ftype: str = None, unit: int = None, extension: str = None +) -> List[Tuple]: + """Get entries from an MF6 namefile. Can select using FTYPE, UNIT, or file extension. + This function only supports MF6 namefiles. + + Parameters + ---------- + path : str + path to a MODFLOW-based model name file + ftype : str + package type + unit : int + file unit number + extension : str + file extension + + Returns + ------- + entries : list of tuples + list of tuples containing FTYPE, UNIT, FNAME, STATUS for each + namefile entry that meets a user-specified value. if no entries + are found, a single-element list containing a None-valued tuple is returned. + + """ + entries = [] + with open(path, "r") as f: + for line in f: + if line.strip() == "": + continue + if line[0] == "#": + continue + ll = line.strip().split() + if len(ll) < 3: + continue + status = "UNKNOWN" + if len(ll) > 3: + status = ll[3].upper() + if ftype is not None: + if ftype.upper() in ll[0].upper(): + filename = os.path.join(os.path.split(path)[0], ll[2]) + entries.append((filename, ll[0], ll[1], status)) + elif unit is not None: + if int(unit) == int(ll[1]): + filename = os.path.join(os.path.split(path)[0], ll[2]) + entries.append((filename, ll[0], ll[1], status)) + elif extension is not None: + filename = os.path.join(os.path.split(path)[0], ll[2]) + ext = os.path.splitext(filename)[1] + if len(ext) > 0: + if ext[0] == ".": + ext = ext[1:] + if extension.lower() == ext.lower(): + entries.append((filename, ll[0], ll[1], status)) + + return entries + + +def get_input_files(namefile): + """Return a list of all the input files in this model. + Parameters + ---------- + namefile : str + path to a MODFLOW-based model name file + Returns + ------- + filelist : list + list of MODFLOW-based model input files + """ + ignore_ext = ( + ".hds", + ".hed", + ".bud", + ".cbb", + ".cbc", + ".ddn", + ".ucn", + ".glo", + ".lst", + ".list", + ".gwv", + ".mv", + ".out", + ) + + srcdir = os.path.dirname(namefile) + filelist = [] + fname = os.path.join(srcdir, namefile) + with open(fname, "r") as f: + lines = f.readlines() + + for line in lines: + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + ext = os.path.splitext(ll[2])[1] + if ext.lower() not in ignore_ext: + if len(ll) > 3: + if "replace" in ll[3].lower(): + continue + filelist.append(ll[2]) + + # Now go through every file and look for other files to copy, + # such as 'OPEN/CLOSE'. If found, then add that file to the + # list of files to copy. + otherfiles = [] + for fname in filelist: + fname = os.path.join(srcdir, fname) + try: + f = open(fname, "r") + for line in f: + + # Skip invalid lines + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + + if "OPEN/CLOSE" in line.upper(): + for i, s in enumerate(ll): + if "OPEN/CLOSE" in s.upper(): + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + otherfiles.append(stmp) + break + except: + print(fname + " does not exist") + + filelist = filelist + otherfiles + + return filelist + + +def get_sim_name(namefiles, rootpth=None): + """Get simulation name. + Parameters + ---------- + namefiles : str or list of strings + path(s) to MODFLOW-based model name files + rootpth : str + optional root directory path (default is None) + Returns + ------- + simname : list + list of namefiles without the file extension + """ + if isinstance(namefiles, str): + namefiles = [namefiles] + sim_name = [] + for namefile in namefiles: + t = namefile.split(os.sep) + if rootpth is None: + idx = -1 + else: + idx = t.index(os.path.split(rootpth)[1]) + + # build dst with everything after the rootpth and before + # the namefile file name. + dst = "" + if idx < len(t): + for d in t[idx + 1 : -1]: + dst += f"{d}_" + + # add namefile basename without extension + dst += t[-1].replace(".nam", "") + sim_name.append(dst) + + return sim_name + + +def get_mf6_nper(tdisfile): + """Return the number of stress periods in the MODFLOW 6 model. + Parameters + ---------- + tdisfile : str + path to the TDIS file + Returns + ------- + nper : int + number of stress periods in the simulation + """ + with open(tdisfile, "r") as f: + lines = f.readlines() + line = [line for line in lines if "NPER" in line.upper()][0] + nper = line.strip().split()[1] + return nper + + +def get_mf6_mshape(disfile): + """Return the shape of the MODFLOW 6 model. + Parameters + ---------- + disfile : str + path to a MODFLOW 6 discretization file + Returns + ------- + mshape : tuple + tuple with the shape of the MODFLOW 6 model. + """ + with open(disfile, "r") as f: + lines = f.readlines() + + d = {} + for line in lines: + + # Skip over blank and commented lines + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + + for key in ["NODES", "NCPL", "NLAY", "NROW", "NCOL"]: + if ll[0].upper() in key: + d[key] = int(ll[1]) + + if "NODES" in d: + mshape = (d["NODES"],) + elif "NCPL" in d: + mshape = (d["NLAY"], d["NCPL"]) + elif "NLAY" in d: + mshape = (d["NLAY"], d["NROW"], d["NCOL"]) + else: + print(d) + raise Exception("Could not determine model shape") + return mshape + + +def get_mf6_files(mfnamefile): + """Return a list of all the MODFLOW 6 input and output files in this model. + Parameters + ---------- + mfnamefile : str + path to the MODFLOW 6 simulation name file + Returns + ------- + filelist : list + list of MODFLOW 6 input files in a simulation + outplist : list + list of MODFLOW 6 output files in a simulation + """ + + srcdir = os.path.dirname(mfnamefile) + filelist = [] + outplist = [] + + filekeys = ["TDIS6", "GWF6", "GWT", "GWF6-GWF6", "GWF-GWT", "IMS6"] + namefilekeys = ["GWF6", "GWT"] + namefiles = [] + + with open(mfnamefile) as f: + + # Read line and skip comments + lines = f.readlines() + + for line in lines: + + # Skip over blank and commented lines + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + + for key in filekeys: + if key in ll[0].upper(): + fname = ll[1] + filelist.append(fname) + + for key in namefilekeys: + if key in ll[0].upper(): + fname = ll[1] + namefiles.append(fname) + + # Go through name files and get files + for namefile in namefiles: + fname = os.path.join(srcdir, namefile) + with open(fname, "r") as f: + lines = f.readlines() + insideblock = False + + for line in lines: + ll = line.upper().strip().split() + if len(ll) < 2: + continue + if ll[0] in "BEGIN" and ll[1] in "PACKAGES": + insideblock = True + continue + if ll[0] in "END" and ll[1] in "PACKAGES": + insideblock = False + + if insideblock: + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + filelist.append(ll[1]) + + # Recursively go through every file and look for other files to copy, + # such as 'OPEN/CLOSE' and 'TIMESERIESFILE'. If found, then + # add that file to the list of files to copy. + flist = filelist + # olist = outplist + while True: + olist = [] + flist, olist = _get_mf6_external_files(srcdir, olist, flist) + # add to filelist + if len(flist) > 0: + filelist = filelist + flist + # add to outplist + if len(olist) > 0: + outplist = outplist + olist + # terminate loop if no additional files + # if len(flist) < 1 and len(olist) < 1: + if len(flist) < 1: + break + + return filelist, outplist + + +def _get_mf6_external_files(srcdir, outplist, files): + """Get list of external files in a MODFLOW 6 simulation. + Parameters + ---------- + srcdir : str + path to a directory containing a MODFLOW 6 simulation + outplist : list + list of output files in a MODFLOW 6 simulation + files : list + list of MODFLOW 6 name files + Returns + ------- + """ + extfiles = [] + + for fname in files: + fname = os.path.join(srcdir, fname) + try: + f = open(fname, "r") + for line in f: + + # Skip invalid lines + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + + if "OPEN/CLOSE" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "OPEN/CLOSE": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + extfiles.append(stmp) + break + + if "TS6" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "FILEIN": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + extfiles.append(stmp) + break + + if "TAS6" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "FILEIN": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + extfiles.append(stmp) + break + + if "OBS6" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "FILEIN": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + extfiles.append(stmp) + break + + if "EXTERNAL" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "EXTERNAL": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + extfiles.append(stmp) + break + + if "FILE" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "FILEIN": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + extfiles.append(stmp) + break + + if "FILE" in line.upper(): + for i, s in enumerate(ll): + if s.upper() == "FILEOUT": + stmp = ll[i + 1] + stmp = stmp.replace('"', "") + stmp = stmp.replace("'", "") + outplist.append(stmp) + break + + except: + print("could not get a list of external mf6 files") + + return extfiles, outplist + + +def get_mf6_ftypes(namefile, ftypekeys): + """Return a list of FTYPES that are in the name file and in ftypekeys. + Parameters + ---------- + namefile : str + path to a MODFLOW 6 name file + ftypekeys : list + list of desired FTYPEs + Returns + ------- + ftypes : list + list of FTYPES that match ftypekeys in namefile + """ + with open(namefile, "r") as f: + lines = f.readlines() + + ftypes = [] + for line in lines: + + # Skip over blank and commented lines + ll = line.strip().split() + if len(ll) < 2: + continue + if line.strip()[0] in ["#", "!"]: + continue + + for key in ftypekeys: + if ll[0].upper() in key: + ftypes.append(ll[0]) + + return ftypes + + +def get_mf6_blockdata(f, blockstr): + """Return list with all non comments between start and end of block + specified by blockstr. + Parameters + ---------- + f : file object + open file object + blockstr : str + name of block to search + Returns + ------- + data : list + list of data in specified block + """ + data = [] + + # find beginning of block + for line in f: + if line[0] != "#": + t = line.split() + if t[0].lower() == "begin" and t[1].lower() == blockstr.lower(): + break + for line in f: + if line[0] != "#": + t = line.split() + if t[0].lower() == "end" and t[1].lower() == blockstr.lower(): + break + else: + data.append(line.rstrip()) + return data