From 86910ed76fcb1796a5c1f25a7c14abfbcc8812fa Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Thu, 22 Jul 2021 14:41:24 -0700 Subject: [PATCH] update(ModflowHob, ModflowDis): Improve HOB file performance (#1158) * 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 --- flopy/modflow/mfdis.py | 70 +++++++++++++++++++++++++++--------------- flopy/modflow/mfhob.py | 18 +++++++++-- 2 files changed, 61 insertions(+), 27 deletions(-) diff --git a/flopy/modflow/mfdis.py b/flopy/modflow/mfdis.py index 7416d26d77..4a47793516 100644 --- a/flopy/modflow/mfdis.py +++ b/flopy/modflow/mfdis.py @@ -289,6 +289,7 @@ def __init__( self.start_datetime = start_datetime # calculate layer thicknesses self.__calculate_thickness() + self._totim = None @property def sr(self): @@ -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): """ @@ -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. @@ -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 ------- @@ -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 @@ -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 @@ -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 ------- @@ -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 diff --git a/flopy/modflow/mfhob.py b/flopy/modflow/mfhob.py index b63c4afc67..1d058e3fc8 100755 --- a/flopy/modflow/mfhob.py +++ b/flopy/modflow/mfhob.py @@ -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() @@ -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 = [] @@ -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( @@ -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