Skip to content

Commit

Permalink
update(ModflowHob, ModflowDis): Improve HOB file performance (#1158)
Browse files Browse the repository at this point in the history 
* reduced repetitive totim array generation for performance improvement
* Added a use_cached_totim flag to get_totim, get_kstp_kper_toffset, get_totim_from_kper_toffset methods in ModflowDis
  • Loading branch information
@jlarsen-usgs
jlarsen-usgs authored Jul 22, 2021
1 parent 78b513a commit 86910ed
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 27 deletions.
70 changes: 46 additions & 24 deletions flopy/modflow/mfdis.py
View file
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
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

0 comments on commit 86910ed

Please sign in to comment.