diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index f420721bd1..e56e20d9c4 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -4,6 +4,7 @@ """ from itertools import repeat +from pprint import pformat import numpy as np import pandas as pd @@ -27,7 +28,7 @@ write_budget, write_head, ) -from flopy.utils.gridutil import uniform_flow_field +from flopy.utils.gridutil import get_disv_kwargs, uniform_flow_field @pytest.fixture @@ -475,47 +476,184 @@ def test_binaryfile_read_context(freyberg_model_path): assert str(e.value) == "seek of closed file", str(e.value) -def test_headfile_reverse_mf6(example_data_path, function_tmpdir): +def test_binaryfile_reverse_mf6_dis(function_tmpdir): + name = "reverse_dis" + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=function_tmpdir, exe_name="mf6" + ) + tdis_rc = [(1, 1, 1.0), (1, 1, 1.0)] + nper = len(tdis_rc) + tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_rc) + ims = flopy.mf6.ModflowIms(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) + dis = flopy.mf6.ModflowGwfdis(gwf, nrow=10, ncol=10) + dis = gwf.get_package("DIS") + nlay = 2 + botm = [1 - (k + 1) for k in range(nlay)] + botm_data = np.array([list(repeat(b, 10 * 10)) for b in botm]).reshape( + (nlay, 10, 10) + ) + dis.nlay = nlay + dis.botm.set_data(botm_data) + ic = flopy.mf6.ModflowGwfic(gwf) + npf = flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True) + chd = flopy.mf6.ModflowGwfchd( + gwf, stress_period_data=[[(0, 0, 0), 1.0], [(0, 9, 9), 0.0]] + ) + budget_file = name + ".bud" + head_file = name + ".hds" + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=budget_file, + head_filerecord=head_file, + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + sim.write_simulation(silent=True) + success, buff = sim.run_simulation(silent=True, report=True) + assert success, pformat(buff) + + # reverse head file in place and check reversal + head_file = flopy.utils.HeadFile(function_tmpdir / head_file, tdis=tdis) + heads = head_file.get_alldata() + assert heads.shape == (nper, 2, 10, 10) + head_file.reverse() + heads_rev = head_file.get_alldata() + assert heads_rev.shape == (nper, 2, 10, 10) + + # reverse budget and write to separate file + budget_file_rev_path = function_tmpdir / f"{budget_file}_rev" + budget_file = flopy.utils.CellBudgetFile( + function_tmpdir / budget_file, tdis=tdis + ) + budget_file.reverse(budget_file_rev_path) + budget_file_rev = flopy.utils.CellBudgetFile( + budget_file_rev_path, tdis=tdis + ) + + for kper in range(nper): + assert np.allclose(heads[kper], heads_rev[-kper + 1]) + budget = budget_file.get_data(text="FLOW-JA-FACE", totim=kper)[0] + budget_rev = budget_file_rev.get_data(text="FLOW-JA-FACE", totim=kper)[ + 0 + ] + assert budget.shape == budget_rev.shape + assert np.allclose(budget, -budget_rev) + + +def test_binaryfile_reverse_mf6_disv(function_tmpdir): + name = "reverse_disv" + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=function_tmpdir, exe_name="mf6" + ) + tdis_rc = [(1, 1, 1.0), (1, 1, 1.0)] + nper = len(tdis_rc) + tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_rc) + ims = flopy.mf6.ModflowIms(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) + dis = flopy.mf6.ModflowGwfdisv( + gwf, **get_disv_kwargs(2, 10, 10, 1.0, 1.0, 25.0, [20.0, 15.0]) + ) + ic = flopy.mf6.ModflowGwfic(gwf) + npf = flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True) + chd = flopy.mf6.ModflowGwfchd( + gwf, stress_period_data=[[(0, 0, 0), 1.0], [(0, 9, 9), 0.0]] + ) + budget_file = name + ".bud" + head_file = name + ".hds" + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=budget_file, + head_filerecord=head_file, + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + sim.write_simulation(silent=True) + success, buff = sim.run_simulation(silent=True) + assert success, pformat(buff) + + # reverse head file in place and check reversal + head_file = flopy.utils.HeadFile(function_tmpdir / head_file, tdis=tdis) + heads = head_file.get_alldata() + assert heads.shape == (nper, 2, 1, 100) + head_file.reverse() + heads_rev = head_file.get_alldata() + assert heads_rev.shape == (nper, 2, 1, 100) + + # reverse budget and write to separate file + budget_file_rev_path = function_tmpdir / f"{budget_file}_rev" + budget_file = flopy.utils.CellBudgetFile( + function_tmpdir / budget_file, tdis=tdis + ) + budget_file.reverse(budget_file_rev_path) + budget_file_rev = flopy.utils.CellBudgetFile( + budget_file_rev_path, tdis=tdis + ) + + for kper in range(nper): + assert np.allclose(heads[kper], heads_rev[-kper + 1]) + budget = budget_file.get_data(text="FLOW-JA-FACE", totim=kper)[0] + budget_rev = budget_file_rev.get_data(text="FLOW-JA-FACE", totim=kper)[ + 0 + ] + assert budget.shape == budget_rev.shape + assert np.allclose(budget, -budget_rev) + + +def test_binaryfile_reverse_mf6_disu(example_data_path, function_tmpdir): # load simulation and extract tdis sim_name = "test006_gwf3" sim = flopy.mf6.MFSimulation.load( sim_name=sim_name, sim_ws=example_data_path / "mf6" / sim_name ) - tdis = sim.get_package("tdis") + tdis_rc = [(1, 1, 1.0), (1, 1, 1.0)] + nper = len(tdis_rc) + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + sim.set_sim_path(function_tmpdir) + sim.write_simulation() + sim.run_simulation() # load head file, providing tdis as kwarg - model_path = example_data_path / "mf6" / sim_name - file_stem = "flow_adj" - file_path = model_path / "expected_output" / f"{file_stem}.hds" - f = HeadFile(file_path, tdis=tdis) - assert isinstance(f, HeadFile) - - # reverse the file - rf_name = f"{file_stem}_rev.hds" - f.reverse(filename=function_tmpdir / rf_name) - rf = HeadFile(function_tmpdir / rf_name) - assert isinstance(rf, HeadFile) + file_path = function_tmpdir / "flow.hds" + head_file = HeadFile(file_path, tdis=tdis) + + # reverse and write to a separate file + head_file_rev_path = function_tmpdir / "flow_rev.hds" + head_file.reverse(filename=head_file_rev_path) + head_file_rev = HeadFile(head_file_rev_path, tdis=tdis) + + # load budget file + file_path = function_tmpdir / "flow.cbc" + budget_file = CellBudgetFile(file_path, tdis=tdis) + + # reverse and write to a separate file + budget_file_rev_path = function_tmpdir / "flow_rev.cbc" + budget_file.reverse(filename=budget_file_rev_path) + budget_file_rev = CellBudgetFile(budget_file_rev_path, tdis=tdis) # check that data from both files have the same shape - assert f.get_alldata().shape == (1, 1, 1, 121) - assert rf.get_alldata().shape == (1, 1, 1, 121) + assert head_file.get_alldata().shape == (nper, 1, 1, 121) + assert head_file_rev.get_alldata().shape == (nper, 1, 1, 121) # check number of records - assert len(f) == 1 - assert len(rf) == 1 + assert len(head_file) == nper + assert len(head_file_rev) == nper + assert len(budget_file) == nper * 2 + assert len(budget_file_rev) == nper * 2 # check that the data are reversed - nrecords = len(f) + nrecords = len(head_file) for idx in range(nrecords - 1, -1, -1): - # check headers - f_header = list(f.recordarray[nrecords - idx - 1]) - rf_header = list(rf.recordarray[idx]) - # todo: these should be equal! + # check headfile headers + f_header = list(head_file.recordarray[nrecords - idx - 1]) + rf_header = list(head_file_rev.recordarray[idx]) assert f_header != rf_header - # check data - f_data = f.get_data(idx=idx)[0] - rf_data = rf.get_data(idx=nrecords - idx - 1)[0] + # check headfile data + f_data = head_file.get_data(idx=idx)[0] + rf_data = head_file_rev.get_data(idx=nrecords - idx - 1)[0] assert f_data.shape == rf_data.shape if f_data.ndim == 1: for row in range(len(f_data)): @@ -525,6 +663,13 @@ def test_headfile_reverse_mf6(example_data_path, function_tmpdir): else: assert np.array_equal(f_data[0][0], rf_data[0][0]) + budget = budget_file.get_data(text="FLOW-JA-FACE", totim=idx)[0] + budget_rev = budget_file_rev.get_data(text="FLOW-JA-FACE", totim=idx)[ + 0 + ] + assert budget.shape == budget_rev.shape + assert np.allclose(budget, -budget_rev) + @pytest.fixture @pytest.mark.mf6 diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index 94074f4275..fb71410926 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile.py @@ -10,6 +10,7 @@ """ import os +import tempfile import warnings from pathlib import Path from typing import List, Optional, Union @@ -20,6 +21,8 @@ from ..utils.datafile import Header, LayerFile from .gridutil import get_lni +HEAD_TEXT = " HEAD" + def write_head( fbin, @@ -28,7 +31,7 @@ def write_head( kper=1, pertim=1.0, totim=1.0, - text=" HEAD", + text=HEAD_TEXT, ilay=1, ): dt = np.dtype( @@ -662,81 +665,99 @@ def __init__( def reverse(self, filename: Optional[os.PathLike] = None): """ - Write a new binary head file with the records in reverse order. - If a new filename is not provided, or if the filename is the same - as the existing filename, the file will be overwritten and data - reloaded from the rewritten/reversed file. + Reverse the time order of the currently loaded binary head file. If a head + file name is not provided or the provided name is the same as the existing + filename, the file will be overwritten and reloaded. Parameters ---------- filename : str or PathLike - Path of the new reversed binary file to create. + Path of the reversed binary head file. """ filename = ( Path(filename).expanduser().absolute() - if filename + if filename is not None else self.filename ) - # header array formats - 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), - ] - ) - - # make sure we have tdis - if self.tdis is None or not any(self.tdis.perioddata.get_data()): - raise ValueError("tdis mu/st be known to reverse head file") - - # extract period data - pd = self.tdis.perioddata.get_data() - - # get maximum period number and total simulation time - kpermx = len(pd) - 1 - tsimtotal = 0.0 - for tpd in pd: - tsimtotal += tpd[0] - - # get total number of records - nrecords = len(self) - - # open backward file - with open(filename, "wb") as fbin: - # loop over head file records in reverse order - for idx in range(nrecords - 1, -1, -1): - # load header array - header = self.recordarray[idx].copy() - - # reverse kstp and kper in the header array - (kstp, kper) = (header["kstp"] - 1, header["kper"] - 1) - kstpmx = pd[kper][1] - 1 - kstpb = kstpmx - kstp - kperb = kpermx - kper - (header["kstp"], header["kper"]) = (kstpb + 1, kperb + 1) + def get_max_kper_kstp_tsim(): + header = self.recordarray[-1] + kper = header["kper"] - 1 + tsim = header["totim"] + kstp = {0: 0} + for i in range(len(self) - 1, -1, -1): + header = self.recordarray[i] + if ( + header["kper"] in kstp + and header["kstp"] > kstp[header["kper"]] + ): + kstp[header["kper"]] += 1 + else: + kstp[header["kper"]] = 0 + return kper, kstp, tsim + + # get max period and time from the head file + maxkper, maxkstp, maxtsim = get_max_kper_kstp_tsim() + # if we have tdis, get max period number and simulation time from it + tdis_maxkper, tdis_maxtsim = None, None + if self.tdis is not None: + pd = self.tdis.perioddata.get_data() + if any(pd): + tdis_maxkper = len(pd) - 1 + tdis_maxtsim = sum([p[0] for p in pd]) + # if we have both, check them against each other + if tdis_maxkper is not None: + assert maxkper == tdis_maxkper, ( + f"Max stress period in binary head file ({maxkper}) != " + f"max stress period in provided tdis ({tdis_maxkper})" + ) + assert maxtsim == tdis_maxtsim, ( + f"Max simulation time in binary head file ({maxtsim}) != " + f"max simulation time in provided tdis ({tdis_maxtsim})" + ) - # reverse totim and pertim in the header array - header["totim"] = tsimtotal - header["totim"] - perlen = pd[kper][0] - header["pertim"] = perlen - header["pertim"] + def reverse_header(header): + """Reverse period, step and time fields in the record header""" + + # reverse kstp and kper headers + kstp = header["kstp"] - 1 + kper = header["kper"] - 1 + header["kstp"] = maxkstp[kper] - kstp + 1 + header["kper"] = maxkper - kper + 1 + + # reverse totim and pertim headers + header["totim"] = maxtsim - header["totim"] + perlen = pd[kper][0] + header["pertim"] = perlen - header["pertim"] + return header + + # reverse record order and write to temporary file + temp_dir_path = Path(tempfile.gettempdir()) + temp_file_path = temp_dir_path / filename.name + with open(temp_file_path, "wb") as f: + for i in range(len(self) - 1, -1, -1): + header = self.recordarray[i].copy() + header = reverse_header(header) + data = self.get_data(idx=i) + ilay = header["ilay"] + write_head( + fbin=f, + data=data[ilay - 1], + kstp=header["kstp"], + kper=header["kper"], + pertim=header["pertim"], + totim=header["totim"], + ilay=ilay, + ) - # write header information - h = np.array(header, dtype=dt) - h.tofile(fbin) + # if we're rewriting the original file, close it first + if filename == self.filename: + self.close() - # load and write data - data = self.get_data(idx=idx)[0][0] - data = np.array(data, dtype=np.float64) - data.tofile(fbin) + # move temp file to destination + temp_file_path.replace(filename) # if we rewrote the original file, reinitialize if filename == self.filename: @@ -2241,21 +2262,23 @@ def close(self): def reverse(self, filename: Optional[os.PathLike] = None): """ - Write a binary cell budget file with the records in reverse order. - If a new filename is not provided, or if the filename is the same - as the existing filename, the file will be overwritten and data - reloaded from the rewritten/reversed file. + Reverse the time order and signs of the currently loaded binary cell budget + file. If a file name is not provided or if the provided name is the same as + the existing filename, the file will be overwritten and reloaded. - Parameters - ---------- + Notes + ----- + While `HeadFile.reverse()` reverses only the temporal order of head data, + this method must reverse not only the order but also the sign (direction) + of the model's intercell flows. filename : str or PathLike, optional - Path of the new reversed binary cell budget file to create. + Path of the reversed binary cell budget file. """ filename = ( Path(filename).expanduser().absolute() - if filename + if filename is not None else self.filename ) @@ -2303,7 +2326,9 @@ def reverse(self, filename: Optional[os.PathLike] = None): nrecords = len(self) # open backward budget file - with open(filename, "wb") as fbin: + temp_dir_path = Path(tempfile.gettempdir()) + temp_file_path = temp_dir_path / filename.name + with open(temp_file_path, "wb") as f: # loop over budget file records in reverse order for idx in range(nrecords - 1, -1, -1): # load header array @@ -2338,7 +2363,7 @@ def reverse(self, filename: Optional[os.PathLike] = None): ] # Note: much of the code below is based on binary_file_writer.py h = np.array(h, dtype=dt1) - h.tofile(fbin) + h.tofile(f) if header["imeth"] == 6: # Write additional header information to the backward budget file h = header[ @@ -2350,7 +2375,7 @@ def reverse(self, filename: Optional[os.PathLike] = None): ] ] h = np.array(h, dtype=dt2) - h.tofile(fbin) + h.tofile(f) # Load data data = self.get_data(idx)[0] data = np.array(data) @@ -2361,7 +2386,7 @@ def reverse(self, filename: Optional[os.PathLike] = None): ndat = len(colnames) - 2 dt = np.dtype([("ndat", np.int32)]) h = np.array([(ndat,)], dtype=dt) - h.tofile(fbin) + h.tofile(f) # Write auxiliary column names naux = ndat - 1 if naux > 0: @@ -2373,12 +2398,12 @@ def reverse(self, filename: Optional[os.PathLike] = None): [(colname, "S16") for colname in colnames[3:]] ) h = np.array(auxtxt, dtype=dt) - h.tofile(fbin) + h.tofile(f) # Write nlist nlist = data.shape[0] dt = np.dtype([("nlist", np.int32)]) h = np.array([(nlist,)], dtype=dt) - h.tofile(fbin) + h.tofile(f) elif header["imeth"] == 1: # Load data data = self.get_data(idx)[0][0][0] @@ -2388,7 +2413,14 @@ def reverse(self, filename: Optional[os.PathLike] = None): else: raise ValueError("not expecting imeth " + header["imeth"]) # Write data - data.tofile(fbin) + data.tofile(f) + + # if we're rewriting the original file, close it first + if filename == self.filename: + self.close() + + # move temp file to destination + temp_file_path.replace(filename) # if we rewrote the original file, reinitialize if filename == self.filename: