diff --git a/modflowapi/extensions/advpaks.py b/modflowapi/extensions/advpaks.py new file mode 100644 index 0000000..f0c6ae0 --- /dev/null +++ b/modflowapi/extensions/advpaks.py @@ -0,0 +1,107 @@ +import numpy as np + +from .data import ListInput +from .pakbase import AdvancedPackage + + +class SfrPakage(AdvancedPackage): + """ + Container for SFR and SFR like packages + + Parameters + ---------- + model : ApiModel + modflowapi model object + pkg_type : str + package type. Ex. "SFR" + pkg_name : str + package name (in the mf6 variables) + sim_package : bool + boolean flag for simulation level packages. Ex. TDIS, IMS + """ + + def __init__(self, model, pkg_type, pkg_name, sim_package=False): + super().__init__(model, pkg_type, pkg_name, sim_package) + + self._diversion_var_arrs = [] + self._set_advanced_variable_addrs("diversions", "_diversion_var_addrs") + self._diversion_vars = ListInput(self, self._diversion_var_arrs, spd=False) + + @property + def diversions(self): + return self._diversion_vars + + @diversions.setter + def diversions(self, recarray): + """ + Setter object to update the diversions data + + """ + if isinstance(recarray, np.recarray): + self._diversion_vars.values = recarray + elif isinstance(recarray, ListInput): + self._diversion_vars.values = recarray.values + elif recarray is None: + self._diversion_vars.values = recarray + else: + raise TypeError(f"{type(recarray)} is not a supported diversions type") + + +class LakPackage(AdvancedPackage): + """ + Container for LAK and LAK like packages + + Parameters + ---------- + model : ApiModel + modflowapi model object + pkg_type : str + package type. Ex. "LAK" + pkg_name : str + package name (in the mf6 variables) + sim_package : bool + boolean flag for simulation level packages. Ex. TDIS, IMS + """ + + def __init__(self, model, pkg_type, pkg_name, sim_package=False): + super().__init__(model, pkg_type, pkg_name, sim_package) + + +class MawPackage(AdvancedPackage): + """ + Container for MAW and MAW like packages + + Parameters + ---------- + model : ApiModel + modflowapi model object + pkg_type : str + package type. Ex. "MAW" + pkg_name : str + package name (in the mf6 variables) + sim_package : bool + boolean flag for simulation level packages. Ex. TDIS, IMS + """ + + def __init__(self, model, pkg_type, pkg_name, sim_package=False): + super().__init__(model, pkg_type, pkg_name, sim_package) + + +class UzfPackage(AdvancedPackage): + """ + Container for UZF and UZF like packages + + Parameters + ---------- + model : ApiModel + modflowapi model object + pkg_type : str + package type. Ex. "UZF" + pkg_name : str + package name (in the mf6 variables) + sim_package : bool + boolean flag for simulation level packages. Ex. TDIS, IMS + """ + + def __init__(self, model, pkg_type, pkg_name, sim_package=False): + super().__init__(model, pkg_type, pkg_name, sim_package) diff --git a/modflowapi/extensions/apimodel.py b/modflowapi/extensions/apimodel.py index 219e3d2..c6c9860 100644 --- a/modflowapi/extensions/apimodel.py +++ b/modflowapi/extensions/apimodel.py @@ -1,9 +1,8 @@ import numpy as np +from .datamodel import get_package_type, gridshape from .pakbase import AdvancedPackage, ArrayPackage, ListPackage, package_factory -gridshape = {"dis": ["nlay", "nrow", "ncol"], "disu": ["nlay", "ncpl"]} - class ApiMbase: """ @@ -15,16 +14,16 @@ class ApiMbase: initialized ModflowApi object name : str modflow model name. ex. "GWF_1", "GWF-GWF_1" - pkg_types : dict - dictionary of package types and ApiPackage class types + pkg_types : None, dict + optional dictionary of package types and ApiPackage class types """ - def __init__(self, mf6, name, pkg_types): + def __init__(self, mf6, name, pkg_types=None): self.mf6 = mf6 self.name = name self._pkg_names = None self._pak_type = None - self.pkg_types = pkg_types + self._pkg_types = pkg_types self.package_dict = {} self._set_package_names() self._create_package_list() @@ -73,10 +72,13 @@ def _create_package_list(self): """ for ix, pkg_name in enumerate(self._pkg_names): pkg_type = self._pak_type[ix].lower() - if pkg_type in self.pkg_types: - basepackage = self.pkg_types[pkg_type] + if self._pkg_types is None: + basepackage = get_package_type(pkg_type) else: - basepackage = AdvancedPackage + if pkg_type in self._pkg_types: + basepackage = self._pkg_types[pkg_type] + else: + basepackage = AdvancedPackage package = package_factory(pkg_type, basepackage) adj_pkg_name = "".join(pkg_type.split("-")) @@ -135,33 +137,6 @@ def __init__(self, mf6, name): else: raise AssertionError(f"Unrecognized discretization type {grid_type}") - pkg_types = { - "dis": ArrayPackage, - "chd": ListPackage, - "drn": ListPackage, - "evt": ListPackage, - "ghb": ListPackage, - "ic": ArrayPackage, - "npf": ArrayPackage, - "rch": ListPackage, - "riv": ListPackage, - "sto": ArrayPackage, - "wel": ListPackage, - # gwt - "dsp": ArrayPackage, - "cnc": ListPackage, - "ist": ArrayPackage, - "mst": ArrayPackage, - "src": ListPackage, - # gwe - "cnd": ArrayPackage, - "est": ArrayPackage, - "cpt": ListPackage, - "esl": ListPackage, - # prt - "mip": ArrayPackage, - } - self.allow_convergence = True self._shape = None self._size = None @@ -169,7 +144,7 @@ def __init__(self, mf6, name): self._usertonode = None self._iteration = 0 - super().__init__(mf6, name, pkg_types) + super().__init__(mf6, name) def __repr__(self): s = f"{self.name}, " diff --git a/modflowapi/extensions/apisimulation.py b/modflowapi/extensions/apisimulation.py index a0c3d01..75b592b 100644 --- a/modflowapi/extensions/apisimulation.py +++ b/modflowapi/extensions/apisimulation.py @@ -311,7 +311,7 @@ def load(mf6): i[:-1].lower() for ix, i in enumerate(mf6.get_value("__INPUT__/SIM/NAM/SLNTYPE")) if idp_names[ix] ] - tmpmdl = ApiMbase(mf6, "", {}) + tmpmdl = ApiMbase(mf6, "") solution_names = list(set(solution_names)) solution_dict = {} for name in solution_names: diff --git a/modflowapi/extensions/data.py b/modflowapi/extensions/data.py index 3e37548..c104da2 100644 --- a/modflowapi/extensions/data.py +++ b/modflowapi/extensions/data.py @@ -3,7 +3,7 @@ import xmipy.errors -class ListInput(object): +class ListInput: """ Data object for storing pointers and working with list based input data @@ -15,40 +15,49 @@ class ListInput(object): optional list of variable addresses mf6 : ModflowApi, None optional ModflowApi object + spd : bool + flag to indicate if the block is loading stress period data or other + list based block data. """ - def __init__(self, parent, var_addrs=None, mf6=None): + def __init__(self, parent, var_addrs=None, mf6=None, spd=True): self.parent = parent + self.var_addrs = var_addrs if self.parent is not None: - self.var_addrs = self.parent.var_addrs + if var_addrs is None: + self.var_addrs = self.parent.var_addrs self.mf6 = self.parent.model.mf6 else: if var_addrs is None or mf6 is None: raise AssertionError("var_addrs and mf6 must be supplied if parent is None") - self.var_addrs = var_addrs + self.mf6 = mf6 self._ptrs = {} + self._compound_ptrs = {} self._nodevars = ("nodelist", "nexg", "maxats") - self._boundvars = ("bound",) + self._bound = "bound" self._maxbound = [0] self._nbound = [0] self._naux = [0] + self._auxvar_name = "auxvar" self._auxnames = [] self._dtype = [] self._reduced_to_var_addr = {} + + self._spd = spd if self.parent._idm_enabled: for var in ("BOUND", "AUXVAR"): self.var_addrs.pop( self.var_addrs.index(self.mf6.get_var_address(var, self.parent.model.name, self.parent.pkg_name)) ) + self.var_addrs.append(self.mf6.get_var_address("AUXVAR_IDM", self.parent.model.name, self.parent.pkg_name)) - self._set_stress_period_data_idm() - else: - self._set_stress_period_data() - def _set_stress_period_data_idm(self): + self._set_data() + + def _set_data(self): """ Method to set stress period data variable pointers to the _ptrs dictionary. Uses IDM updates instead of bound to access variable @@ -56,95 +65,92 @@ def _set_stress_period_data_idm(self): """ # for now we need to add self.parent._bound_vars data to var_addrs for var_addr in self.var_addrs: - try: - values = self.mf6.get_value_ptr(var_addr) - except xmipy.errors.InputError: - if self._naux[0] > 0: - values = self.mf6.get_value(var_addr) - else: + if ":" in var_addr: + addr_pieces = var_addr.split("/") + compound = addr_pieces[-1] + addr_pieces = addr_pieces[:-1] + reduced, ctype, ptr_var = [i.lower() for i in compound.split(":")] + addr_pieces.append(ptr_var.upper()) + var_addr = "/".join(addr_pieces) + + try: + values = self.mf6.get_value_ptr(var_addr) + except xmipy.errors.InputError: continue - reduced = var_addr.split("/")[-1].lower() - if reduced in ("maxbound", "nbound"): - setattr(self, f"_{reduced}", values) - elif reduced in ("nexg", "maxats"): - setattr(self, "_maxbound", values) - setattr(self, "_nbound", values) - elif reduced in ("naux",): - setattr(self, "_naux", values) - elif reduced in ("auxname_cst"): - setattr(self, "_auxnames", list(values)) - else: - self._ptrs[reduced] = values - self._reduced_to_var_addr[reduced] = var_addr - if reduced in self.parent._bound_vars: - typ_str = values.dtype.str - dtype = (reduced, typ_str) - self._dtype.append(dtype) - elif reduced in self._nodevars: - dtype = (reduced, "O") - self._dtype.append(dtype) - elif reduced == "auxvar_idm": - if self._naux == 0: - continue - else: - for ix in range(self._naux[0]): - typ_str = values.dtype.str - dtype = (self._auxnames[ix], typ_str) - self._dtype.append(dtype) + if reduced in ("ndiv",): + nbound = self._special_condition_to_values(ctype, values) + setattr( + self, + "_maxbound", + [ + len(values), + ], + ) + setattr( + self, + "_nbound", + [ + nbound, + ], + ) else: + self._compound_ptrs[ptr_var] = values + self._ptrs[reduced] = f"{ctype}:{ptr_var}" + typ_str = values.dtype.str dtype = (reduced, typ_str) self._dtype.append(dtype) + else: + try: + values = self.mf6.get_value_ptr(var_addr) + except xmipy.errors.InputError: + if self._naux[0] > 0: + values = self.mf6.get_value(var_addr) + else: + continue - def _set_stress_period_data(self): - """ - Method to set stress period data variable pointers to the _ptrs - dictionary - - PENDING DEPRECATION AND REPLACEMENT by _set_stress_period_data_idm() - """ - for var_addr in self.var_addrs: - try: - values = self.mf6.get_value_ptr(var_addr) - except xmipy.errors.InputError: - if self._naux[0] > 0: - values = self.mf6.get_value(var_addr) + reduced = var_addr.split("/")[-1].lower() + if reduced in ("maxbound", "nbound"): + setattr(self, f"_{reduced}", values) + if not self._spd and reduced == "maxbound": + self._nbound = values + elif reduced in ("nexg", "maxats", "nlakes", "nmawwells"): + setattr(self, "_maxbound", values) + setattr(self, "_nbound", values) + elif reduced in ("naux",): + setattr(self, "_naux", values) + elif reduced in ("auxname_cst"): + setattr(self, "_auxnames", list(values)) else: - continue - reduced = var_addr.split("/")[-1].lower() - if reduced in ("maxbound", "nbound"): - setattr(self, f"_{reduced}", values) - elif reduced in ("nexg", "maxats"): - setattr(self, "_maxbound", values) - setattr(self, "_nbound", values) - elif reduced in ("naux",): - setattr(self, "_naux", values) - elif reduced in ("auxname_cst"): - setattr(self, "_auxnames", list(values)) - else: - self._ptrs[reduced] = values - self._reduced_to_var_addr[reduced] = var_addr - if reduced in self._boundvars: - for name in self.parent._bound_vars: + self._ptrs[reduced] = values + self._reduced_to_var_addr[reduced] = var_addr + if reduced == self._bound: + # retain this method for advanced packages + for name in self.parent._bound_vars: + typ_str = values.dtype.str + dtype = (name, typ_str) + self._dtype.append(dtype) + elif reduced in self.parent._bound_vars: typ_str = values.dtype.str - dtype = (name, typ_str) + dtype = (reduced, typ_str) self._dtype.append(dtype) - elif reduced in self._nodevars: - dtype = (reduced, "O") - self._dtype.append(dtype) - elif reduced == "auxvar": - if self._naux == 0: - continue + elif reduced in self._nodevars: + dtype = (reduced, "O") + self._dtype.append(dtype) + elif "auxvar" in reduced: + self._auxvar_name = reduced # == "auxvar_idm": + if self._naux == 0: + continue + else: + for ix in range(self._naux[0]): + typ_str = values.dtype.str + dtype = (self._auxnames[ix], typ_str) + self._dtype.append(dtype) else: - for ix in range(self._naux[0]): - typ_str = values.dtype.str - dtype = (self._auxnames[ix], typ_str) - self._dtype.append(dtype) - else: - typ_str = values.dtype.str - dtype = (reduced, typ_str) - self._dtype.append(dtype) + typ_str = values.dtype.str + dtype = (reduced, typ_str) + self._dtype.append(dtype) def _ptr_to_recarray(self): """ @@ -158,19 +164,29 @@ def _ptr_to_recarray(self): return recarray = np.recarray((self._nbound[0],), self._dtype) for name, ptr in self._ptrs.items(): - if "auxvar" in name and self._naux[0] == 0: + if name == self._auxvar_name and self._naux[0] == 0: continue - values = np.copy(ptr) - if name in self._boundvars: - # note: block slated for deprecation + + if isinstance(ptr, str): + # special condition and not a real ptr + ctype, ptr_name = ptr.split(":") + ptr_vals = self._compound_ptrs[ptr_name] + values = self._special_condition_to_values(ctype, ptr_vals) + else: + values = np.copy(ptr) + + if name == self._bound: + # note: keep block around for advanced packages for ix, nm in enumerate(self.parent._bound_vars): bnd_values = values[0 : self._nbound[0], ix] recarray[nm][0 : self._nbound[0]] = bnd_values + elif name in self.parent._bound_vars and self.parent._idm_enabled: - # new IDM simplification method + # IDM simplification method bnd_values = values[0 : self._nbound[0]].ravel() recarray[name][0 : self._nbound[0]] = bnd_values - elif "auxvar" in name: + + elif name == self._auxvar_name: for ix in range(self._naux[0]): nm = self._auxnames[ix] aux_values = values[0 : self._nbound[0], ix] @@ -186,6 +202,9 @@ def _ptr_to_recarray(self): values = self.parent.model.nodetouser[values] values = list(zip(*np.unravel_index(values, self.parent.model.shape))) + elif name in ("idv", "divreach"): + values -= 1 + values = values[0 : self._nbound[0]] recarray[name][0 : self._nbound[0]] = values @@ -202,11 +221,12 @@ def _recarray_to_ptr(self, recarray): numpy recarray of stress period data """ + if recarray is None: self._nbound[0] = 0 return - if len(recarray) != self._nbound: + if len(recarray) != self._nbound[0]: if len(recarray) > self._maxbound[0]: raise AssertionError( f"Length of stresses ({len(recarray)},) cannot be larger than maxbound value ({self._maxbound[0]},)" @@ -216,33 +236,137 @@ def _recarray_to_ptr(self, recarray): if len(recarray) == 0: return + visited = [] for name in recarray.dtype.names: + if name in visited: + continue if name in self._nodevars: multi_index = tuple(np.array([list(i) for i in recarray[name]]).T) nodes = np.ravel_multi_index(multi_index, self.parent.model.shape) recarray[name] = self.parent.model.usertonode[nodes] + 1 - if name in self.parent._bound_vars: - if "bound" in self._ptrs: - # Block slated for deprecation after IDM inclusion - idx = self.parent._bound_vars.index(name) - bname = "bound" - self._ptrs[bname][0 : self._nbound[0], idx] = recarray[name] - elif self.parent._idm_enabled: - # new IDM simplification - self._ptrs[name][0 : self._nbound[0]] = recarray[name].ravel() - else: - pass + if name in ("divreach",): + self._ptrs[name][0 : self._nbound[0]] = recarray[name].ravel() + 1 + visited.append(name) + + elif self._bound in self._ptrs and name in self.parent._bound_vars: + # Check for bound to support advanced packages + idx = self.parent._bound_vars.index(name) + self._ptrs[self._bound][0 : self._nbound[0], idx] = recarray[name] + visited.append(name) + elif name in self._auxnames: - ptr_name = "auxvar" - if self.parent._idm_enabled: - ptr_name += "_idm" + ptr_name = self._auxvar_name idx = self._auxnames.index(name) self._ptrs[ptr_name][0 : self._nbound[0], idx] = recarray[name] + visited.append(name) + elif name == "auxname_cst": pass + + elif name == self._bound: + pass + else: - self._ptrs[name][0 : self._nbound[0]] = recarray[name].ravel() + if isinstance(self._ptrs[name], str): + visited = self._special_condition_to_ptr(recarray, name, visited) + else: + self._ptrs[name][0 : self._nbound[0]] = recarray[name] + visited.append(name) + + def _special_condition_to_values(self, ctype, inval): + """ + Method to catch and set special, compound conditions where necessary data + is not directly available from MODFLOW + + Parameters + ---------- + ctype : str + condition type + inval : int, float, np.ndarray + input data value + + Returns + ------- + outval : np.array of values + """ + functions = { + "range": lambda x: np.arange(0, int(x), dtype=int), + "count_nonzero": lambda x: np.count_nonzero(x), + "where_idx": lambda x: np.array(np.where(x)[0]), + "where_val": lambda x: x[x != 0], + } + ctype = ctype.lower() + if ctype in ("range",): + inval = inval[0] + + func = functions[ctype] + outval = func(inval) + return outval + + def _special_condition_to_ptr(self, recarray, name, visited): + """ + Method to catch and set special, compound conditions to the associated MODFLOW + ptr where the user data is not directly available from MODFLOW + + Parameters + ---------- + recarray : np.recarray + recarray of user data + name : str + data column name + visited : list + list of visited data columns used to avoid duplicate processing + + + Returns + ------- + visited : list + """ + functions = { + "range": lambda x: [ + len(x), + ], + } + ctype, ptr_name = self._ptrs[name].split(":") + if "where" in ctype: + idx_name = None + val_name = None + if ctype == "where_idx": + idx_name = name + for k, v in self._ptrs.items(): + if isinstance(v, str): + if v == f"where_val:{ptr_name}": + val_name = k + break + + elif ctype == "where_val": + val_name = name + for k, v in self._ptrs.items(): + if isinstance(v, str): + if v == f"where_idx:{ptr_name}": + idx_name = k + break + else: + return + + if idx_name is None or val_name is None: + return + + idx = list(recarray[idx_name].astype(int)) + vals = recarray[val_name].ravel() + if val_name in ("idv",): + vals += 1 + self._compound_ptrs[ptr_name][idx] = vals + visited.extend([idx_name, val_name]) + + else: + func = functions[ctype] + values = func(recarray[name]) + self._compound_ptrs[ptr_name][:] = values[:] + visited.append(name) + + return visited def __getitem__(self, item): recarray = self._ptr_to_recarray() diff --git a/modflowapi/extensions/datamodel.py b/modflowapi/extensions/datamodel.py new file mode 100644 index 0000000..669a0d5 --- /dev/null +++ b/modflowapi/extensions/datamodel.py @@ -0,0 +1,321 @@ +""" +Centralized location to store the "data model"/relationship trees for packages +blocks, and input variables that are used by the modflowapi.extensions code +""" + +gridshape = { + "dis": ["nlay", "nrow", "ncol"], + "disu": [ + "nlay", + "ncpl", + ], +} + + +# Note: HFB variables are not accessible in the memory manager 10/7/2022 +pkgvars = { + "dis": ["top", "bot", "area", "idomain"], + "chd": [ + "nbound", + "maxbound", + "nodelist", + ("bound", ("head",)), + "naux", + "auxname_cst", + "auxvar", + ], + "drn": [ + "nbound", + "maxbound", + "nodelist", + ( + "bound", + ( + "elev", + "cond", + ), + ), + "naux", + "auxname_cst", + "auxvar", + ], + "evt": [ + "nbound", + "maxbound", + "nodelist", + ( + "bound", + ( + "surface", + "rate", + "depth", + ), + ), + # "pxdp:NSEG", "petm:NSEG" + "naux", + "auxname_cst", + "auxvar", + ], + "ghb": [ + "nbound", + "maxbound", + "nodelist", + ( + "bound", + ( + "bhead", + "cond", + ), + ), + "naux", + "auxname_cst", + "auxvar", + ], + "ic": ["strt"], + "npf": ["k11", "k22", "k33", "angle1", "angle2", "angle3", "icelltype"], + "rch": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("recharge",)), + "naux", + "auxname_cst", + "auxvar", + ], + "riv": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("stage", "cond", "rbot")), + "naux", + "auxname_cst", + "auxvar", + ], + "sto": ["iconvert", "ss", "sy"], + "wel": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("q",)), + "naux", + "auxname_cst", + "auxvar", + ], + # gwe model + "cnd": ["alh", "alv", "ath1", "ath2", "atv", "kts"], + "est": ["porosity", "decay", "cps", "rhos"], + "cpt": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("temp",)), + "naux", + "auxname_cst", + "auxvar", + ], + "esl": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("senerrate",)), + "naux", + "auxname_cst", + "auxvar", + ], + # gwt model + "dsp": ["diffc", "alh", "alv", "ath1", "ath2", "atv"], + "cnc": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("conc",)), + "naux", + "auxname_cst", + "auxvar", + ], + "ist": [ + "cim", + "thtaim", + "zetaim", + "decay", + "decay_sorbed", + "bulk_density", + "distcoef", + ], + "mst": ["porosity", "decay", "decay_sorbed", "bulk_density", "distcoef"], + "src": [ + "maxbound", + "nbound", + "nodelist", + ("bound", ("smassrate",)), + "naux", + "auxname_cst", + "auxvar", + ], + # prt model + "mip": ["porosity", "retfactor", "izone"], + # exchange model + "gwf-gwf": ["nexg", "nodem1", "nodem2", "cl1", "cl2", "ihc", "hwva"], + "gwt-gwt": ["nexg", "nodem1", "nodem2", "cl1", "cl2", "ihc", "hwva"], + "gwe-gwe": ["nexg", "nodem1", "nodem2", "cl1", "cl2", "ihc", "hwva"], + # simulation + "ats": [ + "maxats", + "iperats", + "dt0", + "dtmin", + "dtmax", + "dtadj", + "dtfailadj", + ], + "tdis": [ + "nper", + "itmuni", + "kper", + "kstp", + "delt", + "pertim", + "totim,", + "perlen", + "nstp", + "tsmult", + ], + # solution package + "sln-ims": [ + "mxiter", + "dvclose", + "gamma", + "theta", + "akappa", + "amomentum", + "numtrack", + "btol", + "breduc", + "res_lim", + ], + "ims": [ + "niterc", + "dvclose", + "rclose", + "relax", + "ipc", + "droptol", + "north", + "iscl", + "iord", + ], + "sln-ems": [ + "icnvg", + "ttsoln", + ], +} + + +adv_pkgvars = { + "sfr": { + "packagedata": [ + "maxbound", + ( + "ifno:range:maxbound", + "nodelist", + "length", + "width", + "slope", + "strtop", + "bthick", + "hk", + "rough", + "nconnreach", + "ustrf", + "ndiv", + ), + ], + "diversions": [ + "ndiv:count_nonzero:ndiv", + ( + "ifno:where_idx:ndiv", + "idv:where_val:ndiv", + "divreach", # iconr + ), + ], + "perioddata": [ + "maxbound", + "nbound", + ( + "bound", + ("ifno", "sfrsetting", "setting_0", "setting_1"), + ), + ], + }, + "uzf": { + "packagedata": [ + "maxbound", + ( + "ifno:range:maxbound", + "nodelist", + "landflag", + "ivertcon", + "surfdep", + "vks", + "thtr", + "thts", + "thti", + "eps", + ), + ], + "perioddata": [ + "maxbound", + "nbound", + ("bound", ("ifno:range:maxbound", "finf", "pet", "extdp", "extwc", "ha", "hroot", "rootact")), + ], + }, + "lak": { + "packagedata": ["nlakes", ("ifno:range:nlakes", "strt", "nlakeconn")], + }, + "maw": {"packagedata": ["nmawwells", ("ifno:range:nmawwells", "radius", "bot", "strt", "ngwfnodes")]}, +} + + +def get_package_type(pkg_type): + from .advpaks import LakPackage, MawPackage, SfrPakage, UzfPackage + from .pakbase import AdvancedPackage, ArrayPackage, ListPackage, ScalarPackage + + pkg_types = { + "dis": ArrayPackage, + "chd": ListPackage, + "drn": ListPackage, + "evt": ListPackage, + "ghb": ListPackage, + "ic": ArrayPackage, + "npf": ArrayPackage, + "rch": ListPackage, + "riv": ListPackage, + "sto": ArrayPackage, + "wel": ListPackage, + # advanced + "sfr": SfrPakage, + "uzf": UzfPackage, + "lak": LakPackage, + "maw": MawPackage, + # "csub": None, + # gwt + "dsp": ArrayPackage, + "cnc": ListPackage, + "ist": ArrayPackage, + "mst": ArrayPackage, + "src": ListPackage, + # gwe + "cnd": ArrayPackage, + "est": ArrayPackage, + "cpt": ListPackage, + "esl": ListPackage, + # prt + "mip": ArrayPackage, + # sim_level pkgs + "tdis": ScalarPackage, + "ats": ListPackage, + } + if pkg_type in pkg_types: + return pkg_types[pkg_type] + else: + return AdvancedPackage diff --git a/modflowapi/extensions/pakbase.py b/modflowapi/extensions/pakbase.py index fef425f..2a48c9b 100644 --- a/modflowapi/extensions/pakbase.py +++ b/modflowapi/extensions/pakbase.py @@ -1,174 +1,7 @@ import numpy as np from .data import AdvancedInput, ArrayInput, ListInput, ScalarInput - -# Note: HFB variables are not accessible in the memory manager 10/7/2022 -pkgvars = { - "dis": ["top", "bot", "area", "idomain"], - "chd": [ - "nbound", - "maxbound", - "nodelist", - ("bound", ("head",)), - "naux", - "auxname_cst", - "auxvar", - ], - "drn": [ - "nbound", - "maxbound", - "nodelist", - ("bound", ("elev", "cond")), - "naux", - "auxname_cst", - "auxvar", - ], - "evt": [ - "nbound", - "maxbound", - "nodelist", - ("bound", ("surface", "rate", "depth")), - # "pxdp:NSEG", "petm:NSEG" - "naux", - "auxname_cst", - "auxvar", - ], - "ghb": [ - "nbound", - "maxbound", - "nodelist", - ("bound", ("bhead", "cond")), - "naux", - "auxname_cst", - "auxvar", - ], - "ic": ["strt"], - "npf": ["k11", "k22", "k33", "angle1", "angle2", "angle3", "icelltype"], - "rch": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("recharge",)), - "naux", - "auxname_cst", - "auxvar", - ], - "riv": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("stage", "cond", "rbot")), - "naux", - "auxname_cst", - "auxvar", - ], - "sto": ["iconvert", "ss", "sy"], - "wel": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("q",)), - "naux", - "auxname_cst", - "auxvar", - ], - # gwe model - "cnd": ["alh", "alv", "ath1", "ath2", "atv", "kts"], - "est": ["porosity", "decay", "cps", "rhos"], - "cpt": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("temp",)), - "naux", - "auxname_cst", - "auxvar", - ], - "esl": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("senerrate",)), - "naux", - "auxname_cst", - "auxvar", - ], - # gwt model - "dsp": ["diffc", "alh", "alv", "ath1", "ath2", "atv"], - "cnc": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("conc",)), - "naux", - "auxname_cst", - "auxvar", - ], - "ist": [ - "cim", - "thtaim", - "zetaim", - "decay", - "decay_sorbed", - "bulk_density", - "distcoef", - ], - "mst": ["porosity", "decay", "decay_sorbed", "bulk_density", "distcoef"], - "src": [ - "maxbound", - "nbound", - "nodelist", - ("bound", ("smassrate",)), - "naux", - "auxname_cst", - "auxvar", - ], - # prt model - "mip": ["porosity", "retfactor", "izone"], - # exchange model - "gwf-gwf": ["nexg", "nodem1", "nodem2", "cl1", "cl2", "ihc", "hwva"], - "gwt-gwt": ["nexg", "nodem1", "nodem2", "cl1", "cl2", "ihc", "hwva"], - "gwe-gwe": ["nexg", "nodem1", "nodem2", "cl1", "cl2", "ihc", "hwva"], - # simulation - "ats": ["maxats", "iperats", "dt0", "dtmin", "dtmax", "dtadj", "dtfailadj"], - "tdis": [ - "nper", - "itmuni", - "kper", - "kstp", - "delt", - "pertim", - "totim,", - "perlen", - "nstp", - "tsmult", - ], - # solution package - "sln-ims": [ - "mxiter", - "dvclose", - "gamma", - "theta", - "akappa", - "amomentum", - "numtrack", - "btol", - "breduc", - "res_lim", - ], - "ims": [ - "niterc", - "dvclose", - "rclose", - "relax", - "ipc", - "droptol", - "north", - "iscl", - "iord", - ], - "sln-ems": ["icnvg", "ttsoln"], -} +from .datamodel import adv_pkgvars, pkgvars class PackageBase: @@ -232,6 +65,7 @@ def __init__(self, model, pkg_type, pkg_name, child_type, sim_package): for var in self._bound_vars: addr_chk = self.model.mf6.get_var_address(var.upper(), self.model.name, self.pkg_name) if addr_chk in self.model.mf6.get_input_var_names(): + # change this to use idm self._idm_enabled = True var_addrs.append(addr_chk) @@ -629,12 +463,153 @@ class AdvancedPackage(PackageBase): def __init__(self, model, pkg_type, pkg_name, sim_package=False): super().__init__(model, pkg_type, pkg_name.upper(), "advanced", sim_package) + self._idm_enabled = False + self._package_var_addrs = [] + self._sp_var_addrs = [] + self._package_vars = None + self._sp_vars = None + + if pkg_type in adv_pkgvars: + self._adv_var_dict = adv_pkgvars[pkg_type] + + self._set_advanced_variable_addrs("packagedata", "_package_var_addrs") + + if "perioddata" in self._adv_var_dict: + # create variable addresses!!!! + for var in self._adv_var_dict["perioddata"]: + if isinstance(var, tuple): + use_bound = True + for v in var[-1]: + if ":" in v: + use_bound = False + + if use_bound: + self._bound_vars = var[-1] + var = var[0] + else: + for v in var[-1]: + if ":" in v: + tmp = v.split(":")[0] + self._bound_vars.append(tmp) + else: + self._bound_vars.append(v) + var_addr = self._get_advanced_variable_addr(v) + self._sp_var_addrs.append(var_addr) + var = None + + if var is not None: + var_addr = self.model.mf6.get_var_address(var.upper(), self.model.name, self.pkg_name) + self._sp_var_addrs.append(var_addr) + + self._package_vars = ListInput(self, self._package_var_addrs, spd=False) + self._sp_vars = ListInput(self, self._sp_var_addrs) + def __repr__(self): s = f"{self.pkg_type.upper()} Package: {self.pkg_name} \n" s += " Advanced Package, variables only accessible through\n" s += " get_advanced_var() and set_advanced_var() methods" return s + def _set_advanced_variable_addrs(self, block, attr): + """ + General method for setting advanced variable block addresses + to their attributes. Method is used to reduce code duplication + + Parameters + ---------- + block : str + data block key + attr : str + attribute name + + Returns + ------- + None + """ + var_addrs = [] + if block in self._adv_var_dict: + for var in self._adv_var_dict[block]: + if not isinstance(var, tuple): + var_addrs.append(self._get_advanced_variable_addr(var)) + else: + for v in var: + var_addrs.append(self._get_advanced_variable_addr(v)) + + setattr(self, attr, var_addrs) + + def _get_advanced_variable_addr(self, var_str): + """ + Method to create variable addresses for advanced packages that can + include non-standard logic and processing instructions + + Parameters + ---------- + var_str : str + + Returns + ------- + var_addr : str + """ + s = f"{self.model.name}/{self.pkg_name}/{var_str.upper()}" + return s + + @property + def packagedata(self): + """ + Returns a BlockInput object of the packagedata + """ + return self._package_vars + + @packagedata.setter + def packagedata(self, recarray): + """ + Setter method to update the packagedata + + Parameters + ---------- + recarray : np.recarray, ListInput, or None + + """ + if self._package_vars is not None: + if isinstance(recarray, np.recarray): + self._package_vars.values = recarray + elif isinstance(recarray, ListInput): + self._package_vars.values = recarray.values + elif recarray is None: + self._package_vars.values = recarray + else: + raise TypeError(f"{type(recarray)} is not a supported stress_period_data type") + + @property + def maxbound(self): + """ + Returns the "maxbound" value for the stress period + """ + if self._sp_vars is not None: + return self._sp_vars._maxbound[0] + + @property + def stress_period_data(self): + """ + Returns a ListInput object of the current stress_period_data + """ + return self._sp_vars + + @stress_period_data.setter + def stress_period_data(self, recarray): + """ + Setter method to update the current stress_period_data + """ + if self._sp_vars is not None: + if isinstance(recarray, np.recarray): + self._sp_vars.values = recarray + elif isinstance(recarray, ListInput): + self._sp_vars.values = recarray.values + elif recarray is None: + self._sp_vars.values = recarray + else: + raise TypeError(f"{type(recarray)} is not a supported stress_period_data type") + class ApiSlnPackage(ScalarPackage): """