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..e91c01f1bb --- /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 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 = calculate_difftol(h1, h2, htol) + else: + diffmax, indices = calculate_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_concs( + 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 = calculate_difftol(u1, u2, ctol) + else: + diffmax, indices = calculate_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 = calculate_difftol(s1, s2, htol) + else: + diffmax, indices = calculate_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 calculate_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 calculate_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( + 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..9e017225d8 100644 --- a/flopy/utils/mfreadnam.py +++ b/flopy/utils/mfreadnam.py @@ -7,6 +7,7 @@ `_. """ +import os from pathlib import Path, PurePosixPath, PureWindowsPath @@ -266,3 +267,60 @@ def attribs_from_namfile_header(namefile): except: print(f" could not parse start in {namefile}") return defaults + + +def get_entries_from_namefile(namefile, ftype=None, unit=None, extension=None): + """Get entries from a namefile. Can select using FTYPE, UNIT, or file + extension. + + Parameters + ---------- + namefile : 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. + + """ + entries = [] + f = open(namefile, "r") + 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() == ll[0].upper(): + filename = os.path.join(os.path.split(namefile)[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(namefile)[0], ll[2]) + entries.append((filename, ll[0], ll[1], status)) + elif extension is not None: + filename = os.path.join(os.path.split(namefile)[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)) + f.close() + if len(entries) < 1: + entries.append((None, None, None, None)) + return entries