Skip to content

update(ModflowHob, ModflowDis): Improve HOB file performance #1158

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.

Sign up for GitHub

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

Merged
jdhughes-dev merged 1 commit into from
Jul 22, 2021
Merged

update(ModflowHob, ModflowDis): Improve HOB file performance #1158

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
70 changes: 46 additions & 24 deletions flopy/modflow/mfdis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -289,6 +289,7 @@ def __init__(
self.start_datetime = start_datetime
# calculate layer thicknesses
self.__calculate_thickness()
self._totim = None

@property
def sr(self):
Expand Down Expand Up @@ -326,35 +327,45 @@ def checklayerthickness(self):
"""
return (self.parent.modelgrid.thick > 0).all()

def get_totim(self):
def get_totim(self, use_cached=False):
"""
Get the totim at the end of each time step

Parameters
----------
use_cached : bool
method to use cached totim values instead of calculating totim
dynamically


Returns
-------
totim: numpy array
numpy array with simulation totim at the end of each time step

"""
totim = []
nstp = self.nstp.array
perlen = self.perlen.array
tsmult = self.tsmult.array
t = 0.0
for kper in range(self.nper):
m = tsmult[kper]
p = float(nstp[kper])
dt = perlen[kper]
if m > 1:
dt *= (m - 1.0) / (m ** p - 1.0)
else:
dt = dt / p
for kstp in range(nstp[kper]):
t += dt
totim.append(t)
if not use_cached or self._totim is None:
totim = []
nstp = self.nstp.array
perlen = self.perlen.array
tsmult = self.tsmult.array
t = 0.0
for kper in range(self.nper):
m = tsmult[kper]
p = float(nstp[kper])
dt = perlen[kper]
if m > 1:
dt *= m
return np.array(totim, dtype=float)
dt *= (m - 1.0) / (m ** p - 1.0)
else:
dt = dt / p
for kstp in range(nstp[kper]):
t += dt
totim.append(t)
if m > 1:
dt *= m
self._totim = np.array(totim, dtype=float)

return self._totim

def get_final_totim(self):
"""
Expand All @@ -368,7 +379,7 @@ def get_final_totim(self):
"""
return self.get_totim()[-1]

def get_kstp_kper_toffset(self, t=0.0):
def get_kstp_kper_toffset(self, t=0.0, use_cached_totim=False):
"""
Get the stress period, time step, and time offset from passed time.

Expand All @@ -377,6 +388,10 @@ def get_kstp_kper_toffset(self, t=0.0):
t : float
totim to return the stress period, time step, and toffset for
based on time discretization data. Default is 0.
use_cached_totim : bool
optional flag to use a cached calculation of totim, vs. dynamically
calculating totim. Setting to True significantly speeds up looped
operations that call this function (default is False).

Returns
-------
Expand All @@ -391,7 +406,7 @@ def get_kstp_kper_toffset(self, t=0.0):

if t < 0.0:
t = 0.0
totim = self.get_totim()
totim = self.get_totim(use_cached_totim)
nstp = self.nstp.array
ipos = 0
t0 = 0.0
Expand All @@ -415,7 +430,9 @@ def get_kstp_kper_toffset(self, t=0.0):
break
return kstp, kper, toffset

def get_totim_from_kper_toffset(self, kper=0, toffset=0.0):
def get_totim_from_kper_toffset(
self, kper=0, toffset=0.0, use_cached_totim=False
):
"""
Get totim from a passed kper and time offset from the beginning
of a stress period
Expand All @@ -426,6 +443,10 @@ def get_totim_from_kper_toffset(self, kper=0, toffset=0.0):
stress period. Default is 0
toffset : float
time offset relative to the beginning of kper
use_cached_totim : bool
optional flag to use a cached calculation of totim, vs. dynamically
calculating totim. Setting to True significantly speeds up looped
operations that call this function (default is False).

Returns
-------
Expand All @@ -443,8 +464,9 @@ def get_totim_from_kper_toffset(self, kper=0, toffset=0.0):
+ "must be less than "
+ "to nper ({}).".format(self.nper)
)
raise ValueError()
totim = self.get_totim()
raise ValueError(msg)

totim = self.get_totim(use_cached_totim)
nstp = self.nstp.array
ipos = 0
t0 = 0.0
Expand Down
18 changes: 15 additions & 3 deletions flopy/modflow/mfhob.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -376,6 +376,10 @@ def load(cls, f, model, ext_unit_dict=None, check=True):

# read datasets 3-6
nobs = 0

# set to False for 1st call to ensure that totim cache is updated
use_cached_totim = False

while True:
# read dataset 3
line = f.readline()
Expand Down Expand Up @@ -429,11 +433,12 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
itt = 1
irefsp0 -= 1
totim = model.dis.get_totim_from_kper_toffset(
irefsp0, toffset * tomulth
irefsp0, toffset * tomulth, use_cached_totim
)
names = [obsnam]
tsd = [totim, hob]
nobs += 1
use_cached_totim = True
else:
names = []
tsd = []
Expand All @@ -449,11 +454,12 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
irefsp = int(t[1]) - 1
toffset = float(t[2])
totim = model.dis.get_totim_from_kper_toffset(
irefsp, toffset * tomulth
irefsp, toffset * tomulth, use_cached_totim
)
hob = float(t[3])
tsd.append([totim, hob])
nobs += 1
use_cached_totim = True

obs_data.append(
HeadObservation(
Expand Down Expand Up @@ -680,16 +686,22 @@ def __init__(
)
raise ValueError(msg)

# set use_cached_totim to False first to ensure totim is updated
use_cached_totim = False

# create time_series_data
self.time_series_data = self._get_empty(ncells=shape[0])
for idx in range(self.nobs):
t = time_series_data[idx, 0]
kstp, kper, toffset = model.dis.get_kstp_kper_toffset(t)
kstp, kper, toffset = model.dis.get_kstp_kper_toffset(
t, use_cached_totim
)
self.time_series_data[idx]["totim"] = t
self.time_series_data[idx]["irefsp"] = kper
self.time_series_data[idx]["toffset"] = toffset / tomulth
self.time_series_data[idx]["hobs"] = time_series_data[idx, 1]
self.time_series_data[idx]["obsname"] = names[idx]
use_cached_totim = True

if self.nobs > 1:
self.irefsp = -self.nobs
Expand Down