From 64634fec80095513d1f340491ae63e8de4d2a1fc Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Thu, 22 Jun 2023 08:16:24 -0400 Subject: [PATCH] refactor(binaryfile): require tdis for reverse() methods --- autotest/test_binaryfile.py | 45 +++++++++++++++---- flopy/utils/binaryfile.py | 86 ++++++++++++++++++++----------------- flopy/utils/datafile.py | 2 + 3 files changed, 85 insertions(+), 48 deletions(-) diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index 12dea639a3..1982b620c6 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -6,6 +6,7 @@ from matplotlib import pyplot as plt from matplotlib.axes import Axes +from flopy.mf6.modflow import MFSimulation from flopy.modflow import Modflow from flopy.utils import ( BinaryHeader, @@ -552,9 +553,17 @@ def test_cellbudgetfile_readrecord_waux(example_data_path): reason="failing, need to modify CellBudgetFile.reverse to support mf2005?" ) def test_cellbudgetfile_reverse_mf2005(example_data_path, function_tmpdir): - mf2005_model_path = example_data_path / "mf2005_test" - cbc_fname = mf2005_model_path / "test1tr.gitcbc" - f = CellBudgetFile(cbc_fname) + sim_name = "test1tr" + + # load simulation and extract tdis + sim = MFSimulation.load( + sim_name=sim_name, sim_ws=example_data_path / "mf2005_test" + ) + tdis = sim.get_package("tdis") + + mf2005_model_path = example_data_path / sim_name + cbc_fname = mf2005_model_path / f"{sim_name}.gitcbc" + f = CellBudgetFile(cbc_fname, tdis=tdis) assert isinstance(f, CellBudgetFile) rf_name = "test1tr_rev.gitcbc" @@ -564,15 +573,24 @@ def test_cellbudgetfile_reverse_mf2005(example_data_path, function_tmpdir): def test_cellbudgetfile_reverse_mf6(example_data_path, function_tmpdir): - model_path = example_data_path / "mf6" / "test006_gwf3" + sim_name = "test006_gwf3" + + # load simulation and extract tdis + sim = MFSimulation.load( + sim_name=sim_name, sim_ws=example_data_path / "mf6" / sim_name + ) + tdis = sim.get_package("tdis") + + # load cell budget 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}.cbc" - f = CellBudgetFile(file_path) + f = CellBudgetFile(file_path, tdis=tdis) assert isinstance(f, CellBudgetFile) # reverse the file rf_name = f"{file_stem}_rev.cbc" - f.reverse(function_tmpdir / rf_name) + f.reverse(filename=function_tmpdir / rf_name) rf = CellBudgetFile(function_tmpdir / rf_name) assert isinstance(rf, CellBudgetFile) @@ -603,15 +621,24 @@ def test_cellbudgetfile_reverse_mf6(example_data_path, function_tmpdir): def test_headfile_reverse_mf6(example_data_path, function_tmpdir): - model_path = example_data_path / "mf6" / "test006_gwf3" + sim_name = "test006_gwf3" + + # load simulation and extract tdis + sim = MFSimulation.load( + sim_name=sim_name, sim_ws=example_data_path / "mf6" / sim_name + ) + tdis = sim.get_package("tdis") + + # load cell budget 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) + f = HeadFile(file_path, tdis=tdis) assert isinstance(f, HeadFile) # reverse the file rf_name = f"{file_stem}_rev.hds" - f.reverse(function_tmpdir / rf_name) + f.reverse(filename=function_tmpdir / rf_name) rf = HeadFile(function_tmpdir / rf_name) assert isinstance(rf, HeadFile) diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index 1d120f118f..a669a8a203 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile.py @@ -11,7 +11,7 @@ import os import warnings from pathlib import Path -from typing import Union +from typing import Optional, Union import numpy as np @@ -636,7 +636,7 @@ def __init__( ) super().__init__(filename, precision, verbose, kwargs) - def reverse(self, filename: os.PathLike): + def reverse(self, filename: Optional[os.PathLike] = None): """ Write a new binary file with the records in reverse order. @@ -661,46 +661,48 @@ def reverse(self, filename: os.PathLike): ] ) + # 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 = max(self.recordarray["kper"]) - tsimtotal = max(self.get_times()) + kpermx = len(pd) - 1 + tsimtotal = 0.0 + for tpd in pd: + tsimtotal += tpd[0] # get total number of records - nrecords = len(self.get_alldata()) - - # reverse kstpkper and times - kstpkper = list(reversed(self.kstpkper.copy())) - times = list(reversed(self.times.copy())) + nrecords = self.recordarray.shape[0] # open backward file - with open(filename, "wb") as fbin: + with open(filename if filename else self.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"], header["kper"]) - kstpmx = next(iter([kk for kk in kstpkper if kk[1] == kper]))[ - 0 - ] - (header["kstp"], header["kper"]) = ( - kstpmx - kstp + 1, - kpermx - kper + 1, - ) + (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) # reverse totim and pertim in the header array header["totim"] = tsimtotal - header["totim"] - perlen = times[kper - 1] - times[kper - 2] + perlen = pd[kper][0] header["pertim"] = perlen - header["pertim"] - # write header information to the backward head file + # write header information h = np.array(header, dtype=dt) h.tofile(fbin) - # load data + + # load and write data data = self.get_data(idx=idx)[0][0] data = np.array(data, dtype=np.float64) - # write data data.tofile(fbin) @@ -845,6 +847,8 @@ def __init__( if "dis" in kwargs.keys(): self.dis = kwargs.pop("dis") self.modelgrid = self.dis.parent.modelgrid + if "tdis" in kwargs.keys(): + self.tdis = kwargs.pop("tdis") if "sr" in kwargs.keys(): from ..discretization import StructuredGrid, UnstructuredGrid @@ -2018,9 +2022,8 @@ def close(self): Close the file handle """ self.file.close() - return - def reverse(self, filename: os.PathLike): + def reverse(self, filename: Optional[os.PathLike] = None): """ Write a new binary cell budget file with the records in reverse order. @@ -2055,17 +2058,25 @@ def reverse(self, filename: os.PathLike): ] ) + # make sure we have tdis + if self.tdis is None or not any(self.tdis.perioddata.get_data()): + raise ValueError( + "tdis must be known to reverse a cell budget file" + ) + + # extract perioddata + pd = self.tdis.perioddata.get_data() + # get maximum period number and total simulation time - kpermx = self.nper - 1 - tsimtotal = max(self.times) + nper = len(pd) + kpermx = nper - 1 + tsimtotal = 0.0 + for tpd in pd: + tsimtotal += tpd[0] # get number of records nrecords = self.get_nrecords() - # reverse kstpkper - kstpkper = list(reversed(self.kstpkper.copy())) - times = list(reversed(self.times.copy())) - # open backward budget file with open(filename, "wb") as fbin: # loop over budget file records in reverse order @@ -2074,18 +2085,15 @@ def reverse(self, filename: os.PathLike): header = self.recordarray[idx] # reverse kstp and kper in the header array - (kstp, kper) = (header["kstp"], header["kper"]) - kstpmx = next(iter([kk for kk in kstpkper if kk[1] == kper]))[ - 0 - ] - (header["kstp"], header["kper"]) = ( - kstpmx - kstp, - kpermx - kper, - ) + (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) # reverse totim and pertim in the header array header["totim"] = tsimtotal - header["totim"] - perlen = times[kper - 1] - times[kper - 2] + perlen = pd[kper][0] header["pertim"] = perlen - header["pertim"] # Write main header information to backward budget file diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 73a7f61b66..cbe0ac4410 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -196,6 +196,8 @@ def __init__( if "dis" in kwargs.keys(): self.dis = kwargs.pop("dis") self.mg = self.dis.parent.modelgrid + if "tdis" in kwargs.keys(): + self.tdis = kwargs.pop("tdis") if "modelgrid" in kwargs.keys(): self.mg = kwargs.pop("modelgrid") if len(kwargs.keys()) > 0: