Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat(binaryfile): add head/budget file reversal script #2383

Merged
merged 3 commits into from
Nov 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions autotest/test_binaryfile.py
Original file line number Diff line number Diff line change
Expand Up @@ -502,7 +502,7 @@ def test_binaryfile_reverse_mf6_dis(function_tmpdir):
assert success, pformat(buff)

# reverse head file in place and check reversal
head_file = flopy.utils.HeadFile(function_tmpdir / head_file, tdis=tdis)
head_file = flopy.utils.HeadFile(function_tmpdir / head_file)
heads = head_file.get_alldata()
assert heads.shape == (nper, 2, 10, 10)
head_file.reverse()
Expand All @@ -511,9 +511,9 @@ def test_binaryfile_reverse_mf6_dis(function_tmpdir):

# 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 = flopy.utils.CellBudgetFile(function_tmpdir / budget_file)
budget_file.reverse(budget_file_rev_path)
budget_file_rev = flopy.utils.CellBudgetFile(budget_file_rev_path, tdis=tdis)
budget_file_rev = flopy.utils.CellBudgetFile(budget_file_rev_path)

for kper in range(nper):
assert np.allclose(heads[kper], heads_rev[-kper + 1])
Expand Down Expand Up @@ -554,7 +554,7 @@ def test_binaryfile_reverse_mf6_disv(function_tmpdir):
assert success, pformat(buff)

# reverse head file in place and check reversal
head_file = flopy.utils.HeadFile(function_tmpdir / head_file, tdis=tdis)
head_file = flopy.utils.HeadFile(function_tmpdir / head_file)
heads = head_file.get_alldata()
assert heads.shape == (nper, 2, 1, 100)
head_file.reverse()
Expand All @@ -563,7 +563,7 @@ def test_binaryfile_reverse_mf6_disv(function_tmpdir):

# 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 = flopy.utils.CellBudgetFile(function_tmpdir / budget_file)
budget_file.reverse(budget_file_rev_path)
budget_file_rev = flopy.utils.CellBudgetFile(budget_file_rev_path, tdis=tdis)

Expand Down Expand Up @@ -595,7 +595,7 @@ def test_binaryfile_reverse_mf6_disu(example_data_path, function_tmpdir):
# 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)
head_file_rev = HeadFile(head_file_rev_path)

# load budget file
file_path = function_tmpdir / "flow.cbc"
Expand All @@ -604,7 +604,7 @@ def test_binaryfile_reverse_mf6_disu(example_data_path, function_tmpdir):
# 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)
budget_file_rev = CellBudgetFile(budget_file_rev_path)

# check that data from both files have the same shape
assert head_file.get_alldata().shape == (nper, 1, 1, 121)
Expand Down
96 changes: 48 additions & 48 deletions flopy/utils/binaryfile.py → flopy/utils/binaryfile/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@
import numpy as np
import pandas as pd

from ..utils.datafile import Header, LayerFile
from .gridutil import get_lni
from ..datafile import Header, LayerFile
from ..gridutil import get_lni

HEAD_TEXT = " HEAD"

Expand Down Expand Up @@ -663,38 +663,28 @@ def get_max_kper_kstp_tsim():
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})"
)
prev_kper = None
perlen = None

def reverse_header(header):
"""Reverse period, step and time fields in the record header"""

nonlocal prev_kper
nonlocal perlen

# reverse kstp and kper headers
kstp = header["kstp"] - 1
kper = header["kper"] - 1
header["kstp"] = maxkstp[kper] - kstp + 1
header["kper"] = maxkper - kper + 1

if kper != prev_kper:
perlen = header["pertim"]
prev_kper = kper

# reverse totim and pertim headers
header["totim"] = maxtsim - header["totim"]
perlen = pd[kper][0]
header["pertim"] = perlen - header["pertim"]
return header

Expand Down Expand Up @@ -1022,7 +1012,6 @@ def __init__(
self.paknamlist_from = []
self.paknamlist_to = []
self.compact = True # compact budget file flag

self.dis = None
self.modelgrid = None
if "model" in kwargs.keys():
Expand Down Expand Up @@ -2237,24 +2226,46 @@ def reverse(self, filename: Optional[os.PathLike] = None):
]
)

# 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")
nrecords = len(self)
target = filename

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

maxkper, maxkstp, maxtsim = get_max_kper_kstp_tsim()
prev_kper = None
perlen = None

# extract perioddata
pd = self.tdis.perioddata.get_data()
def reverse_header(header):
"""Reverse period, step and time fields in the record header"""

# get maximum period number and total simulation time
nper = len(pd)
kpermx = nper - 1
tsimtotal = 0.0
for tpd in pd:
tsimtotal += tpd[0]
nonlocal prev_kper
nonlocal perlen

# get number of records
nrecords = len(self)
# reverse kstp and kper headers
kstp = header["kstp"] - 1
kper = header["kper"] - 1
header["kstp"] = maxkstp[kper] - kstp + 1
header["kper"] = maxkper - kper + 1

target = filename
if kper != prev_kper:
perlen = header["pertim"] - 1
prev_kper = kper

# reverse totim and pertim headers
header["totim"] = maxtsim - header["totim"]
header["pertim"] = perlen - header["pertim"]
return header

# if rewriting the same file, write
# temp file then copy it into place
Expand All @@ -2269,18 +2280,7 @@ def reverse(self, filename: Optional[os.PathLike] = None):
for idx in range(nrecords - 1, -1, -1):
# load header array
header = self.recordarray[idx]

# 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)

# reverse totim and pertim in the header array
header["totim"] = tsimtotal - header["totim"]
perlen = pd[kper][0]
header["pertim"] = perlen - header["pertim"]
header = reverse_header(header)

# Write main header information to backward budget file
h = header[
Expand Down
31 changes: 31 additions & 0 deletions flopy/utils/binaryfile/reverse.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
import argparse
from pathlib import Path

from flopy.utils.binaryfile import CellBudgetFile, HeadFile

if __name__ == "__main__":
"""Reverse head or budget files."""

parser = argparse.ArgumentParser(description="Reverse head or budget files.")
parser.add_argument(
"--infile",
"-i",
type=str,
help="Input file.",
)
parser.add_argument(
"--outfile",
"-o",
type=str,
help="Output file.",
)
args = parser.parse_args()
infile = Path(args.infile)
outfile = Path(args.outfile)
suffix = infile.suffix.lower()
if suffix in [".hds", ".hed"]:
HeadFile(infile).reverse(outfile)
elif suffix in [".bud", ".cbc"]:
CellBudgetFile(infile).reverse(outfile)
else:
raise ValueError(f"Unrecognized file suffix: {suffix}")
Loading