Skip to content

Commit

Permalink
feat(binaryfile): add head/budget file reversal script
Browse files Browse the repository at this point in the history
  • Loading branch information
wpbonelli committed Nov 26, 2024
1 parent 125b926 commit eef5d33
Show file tree
Hide file tree
Showing 4 changed files with 80 additions and 49 deletions.
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"]
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}")
2 changes: 1 addition & 1 deletion flopy/utils/util_array.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,8 @@
import numpy as np

from ..datbase import DataInterface, DataType
from ..utils.binaryfile import BinaryHeader
from ..utils.flopy_io import line_parse
from .binaryfile import BinaryHeader


class ArrayFormat:
Expand Down
Binary file added freyberg-rev.hds
Binary file not shown.

0 comments on commit eef5d33

Please sign in to comment.