From 7957756f08b54a838e9d99865ad1e61584831224 Mon Sep 17 00:00:00 2001 From: Hua Zhang Date: Fri, 22 Nov 2024 09:47:59 -0500 Subject: [PATCH] Feat (mfusg): FloPy MODFLOW-USG support (with .bas file example) #1476 * Feat(mfusg): add mfusg transport version 2.4.0 support fixes and updates included with feature: New Features * MfUsgBct : Block Centered Transport (BCT) * MfUsgPcb : Prescribed Concentration Boundary (PCB) * MfUsgDdf: Density Driven Flow (DDF) * MfUsgDpf: Dual Porosity Flow (DPF) * MfUsgDpt: Dual Porosity Transport (DPT) * MfUsgMdt: Matrix Diffusion Transport (MDT) * MfUsgRch: Recharge (RCH) * MfUsgLak: Lake (LAK) Modified Packages * MfUsgBcf : Block Centered Flow (BCF) * MfUsgCln : Connected Linear Network (CLN) * MfUsgSms : Sparse Matrix Solver (SMS) * MfUsgWel : Well (WEL) * ModflowFhb : Flow and Head Boundary (FHB) See Panday, S., 2024; USG-Transport Version 2.4.0: Transport and Other Enhancements to MODFLOW-USG, GSI Environmental, July 2024 http://www.gsi-net.com/en/software/free-software/USG-Transport.html --- flopy/mfusg/__init__.py | 18 + flopy/mfusg/mfusg.py | 19 +- flopy/mfusg/mfusgbcf.py | 62 +- flopy/mfusg/mfusgbct.py | 1092 +++++++++++++++++++++++++++++++++ flopy/mfusg/mfusgcln.py | 185 +++++- flopy/mfusg/mfusgddf.py | 219 +++++++ flopy/mfusg/mfusgdisu.py | 2 + flopy/mfusg/mfusgdpf.py | 327 ++++++++++ flopy/mfusg/mfusgdpt.py | 616 +++++++++++++++++++ flopy/mfusg/mfusglak.py | 924 ++++++++++++++++++++++++++++ flopy/mfusg/mfusgmdt.py | 454 ++++++++++++++ flopy/mfusg/mfusgoc.py | 1240 ++++++++++++++++++++++++++++++++++++++ flopy/mfusg/mfusgpcb.py | 294 +++++++++ flopy/mfusg/mfusgrch.py | 518 ++++++++++++++++ flopy/mfusg/mfusgsms.py | 108 ++-- flopy/mfusg/mfusgwel.py | 4 - flopy/modflow/mffhb.py | 5 +- 17 files changed, 6003 insertions(+), 84 deletions(-) create mode 100644 flopy/mfusg/mfusgbct.py create mode 100644 flopy/mfusg/mfusgddf.py create mode 100644 flopy/mfusg/mfusgdpf.py create mode 100644 flopy/mfusg/mfusgdpt.py create mode 100644 flopy/mfusg/mfusglak.py create mode 100644 flopy/mfusg/mfusgmdt.py create mode 100644 flopy/mfusg/mfusgoc.py create mode 100644 flopy/mfusg/mfusgpcb.py create mode 100644 flopy/mfusg/mfusgrch.py diff --git a/flopy/mfusg/__init__.py b/flopy/mfusg/__init__.py index 6790a64896..f657fd3c56 100644 --- a/flopy/mfusg/__init__.py +++ b/flopy/mfusg/__init__.py @@ -3,10 +3,19 @@ from .cln_dtypes import MfUsgClnDtypes from .mfusg import MfUsg from .mfusgbcf import MfUsgBcf +from .mfusgbct import MfUsgBct from .mfusgcln import MfUsgCln +from .mfusgddf import MfUsgDdf from .mfusgdisu import MfUsgDisU +from .mfusgdpf import MfUsgDpf +from .mfusgdpt import MfUsgDpt from .mfusggnc import MfUsgGnc +from .mfusglak import MfUsgLak from .mfusglpf import MfUsgLpf +from .mfusgmdt import MfUsgMdt +from .mfusgoc import MfUsgOc +from .mfusgpcb import MfUsgPcb +from .mfusgrch import MfUsgRch from .mfusgsms import MfUsgSms from .mfusgwel import MfUsgWel @@ -20,4 +29,13 @@ "MfUsgClnDtypes", "MfUsgSms", "MfUsgGnc", + "MfUsgBct", + "MfUsgPcb", + "MfUsgDdf", + "MfUsgMdt", + "MfUsgDpf", + "MfUsgDpt", + "MfUsgRch", + "MfUsgOc", + "MfUsgLak", ] diff --git a/flopy/mfusg/mfusg.py b/flopy/mfusg/mfusg.py index 5dc926483f..e82dade8be 100644 --- a/flopy/mfusg/mfusg.py +++ b/flopy/mfusg/mfusg.py @@ -91,6 +91,11 @@ def __init__( verbose=self.verbose, **kwargs, ) + + self.mcomp = 0 # number of chemical components + self.iheat = 0 # flag for heat transport + self.idpf = 0 # flag for dual porosity flow + # Create a dictionary to map package with package object. # This is used for loading models. self.mfnam_packages = { @@ -104,15 +109,12 @@ def __init__( "fhb": flopy.modflow.ModflowFhb, "drn": flopy.modflow.ModflowDrn, "drt": flopy.modflow.ModflowDrt, - "rch": flopy.modflow.ModflowRch, "evt": flopy.modflow.ModflowEvt, "ghb": flopy.modflow.ModflowGhb, "riv": flopy.modflow.ModflowRiv, "str": flopy.modflow.ModflowStr, "sfr": flopy.modflow.ModflowSfr2, - "lak": flopy.modflow.ModflowLak, "gage": flopy.modflow.ModflowGage, - "oc": flopy.modflow.ModflowOc, "sub": flopy.modflow.ModflowSub, "swt": flopy.modflow.ModflowSwt, "disu": flopy.mfusg.MfUsgDisU, @@ -122,6 +124,15 @@ def __init__( "lpf": flopy.mfusg.MfUsgLpf, "cln": flopy.mfusg.MfUsgCln, "gnc": flopy.mfusg.MfUsgGnc, + "bct": flopy.mfusg.MfUsgBct, + "pcb": flopy.mfusg.MfUsgPcb, + "ddf": flopy.mfusg.MfUsgDdf, + "mdt": flopy.mfusg.MfUsgMdt, + "dpf": flopy.mfusg.MfUsgDpf, + "dpt": flopy.mfusg.MfUsgDpt, + "rch": flopy.mfusg.MfUsgRch, + "oc": flopy.mfusg.MfUsgOc, + "lak": flopy.mfusg.MfUsgLak, } def __repr__(self): @@ -529,7 +540,7 @@ def fmt_string(array): if vtype in ("i", "b"): fmts.append("%10d") elif vtype == "f": - fmts.append("%14.6g") + fmts.append("%10.2g") elif vtype == "o": fmts.append("%10s") elif vtype == "s": diff --git a/flopy/mfusg/mfusgbcf.py b/flopy/mfusg/mfusgbcf.py index ae0edc2b80..8200afbb1f 100644 --- a/flopy/mfusg/mfusgbcf.py +++ b/flopy/mfusg/mfusgbcf.py @@ -138,6 +138,12 @@ def __init__( ihdwet=0, ikvflag=0, ikcflag=0, + tabrich=False, + nuzones=0, + nutabrows=0, + bubblept=False, + fullydry=False, + altsto=False, tran=1.0, hy=1.0, vcont=1.0, @@ -194,6 +200,13 @@ def __init__( self.ikvflag = ikvflag self.ikcflag = ikcflag + self.tabrich = tabrich + self.nuzones = nuzones + self.nutabrows = nutabrows + self.bubblept = bubblept + self.fullydry = fullydry + self.altsto = altsto + self.kv = kv self.anglex = anglex self.ksat = ksat @@ -255,9 +268,19 @@ def write_file(self, f=None): f_obj.write( f" {self.ipakcb:9d} {self.hdry:9.3G} {self.iwdflg:9d}" f" {self.wetfct:9.3G} {self.iwetit:9d} {self.ihdwet:9d}" - f" {self.ikvflag:9d} {self.ikcflag:9d}\n" + f" {self.ikvflag:9d} {self.ikcflag:9d}" ) + if self.tabrich: + f_obj.write(f" TABRICH {self.nuzones:9d} {self.nutabrows:9d}") + if self.bubblept: + f_obj.write(" BUBBLEPT") + if self.fullydry: + f_obj.write(" FULLYDRY") + if self.altsto: + f_obj.write(" ALTSTO") + f_obj.write("\n") + # LAYCON array for layer in range(nlay): if self.intercellt[layer] > 0: @@ -386,6 +409,37 @@ def load(cls, f, model, ext_unit_dict=None): ikvflag = type_from_iterable(text_list, index=6, _type=int, default_val=0) ikcflag = type_from_iterable(text_list, index=7, _type=int, default_val=0) + # options + if "TABRICH" in text_list: + idx = text_list.index("TABRICH") + tabrich = True + nuzones = float(text_list[idx + 1]) + nutabrows = float(text_list[idx + 2]) + else: + tabrich = False + nuzones = None + nutabrows = None + + if "BUBBLEPT" in text_list: + bubblept = True + else: + bubblept = False + + if "FULLYDRY" in text_list: + fullydry = True + else: + fullydry = False + + if "ALTSTO" in text_list: + altsto = True + else: + altsto = False + + # item 1c and d -- Not implemented + if tabrich: + iuzontab = [] + retcrvs = [] + # LAYCON array laycon, intercellt = cls._load_laycon(f_obj, model) @@ -446,6 +500,12 @@ def load(cls, f, model, ext_unit_dict=None): ihdwet=ihdwet, ikvflag=ikvflag, ikcflag=ikcflag, + tabrich=tabrich, + nuzones=nuzones, + nutabrows=nutabrows, + bubblept=bubblept, + fullydry=fullydry, + altsto=altsto, tran=tran, hy=hy, vcont=vcont, diff --git a/flopy/mfusg/mfusgbct.py b/flopy/mfusg/mfusgbct.py new file mode 100644 index 0000000000..96636ff0a7 --- /dev/null +++ b/flopy/mfusg/mfusgbct.py @@ -0,0 +1,1092 @@ +""" +Mfusgbct module. + +Contains the MfUsgBct class. Note that the user can +access the MfUsgBct class as `flopy.mfusg.MfUsgBct`. +""" + +import numpy as np + +from ..pakbase import Package +from ..utils import Util2d, Util3d, read1d +from ..utils.utils_def import ( + get_open_file_object, + get_util2d_shape_for_layer, + type_from_iterable, +) +from .mfusg import MfUsg + + +class MfUsgBct(Package): + """Block Centered Transport (BCT) Package Class for MODFLOW-USG Transport. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.modflow.Modflow`) to which + this package will be added. + itrnsp : int (0,1,2,3,4,5,-1), (default is 1) + transport simulation flag + ipakcb : int (0,1,-1), (default is 0) + a flag and a unit number >0 . + mcomp : int (default is 1) + number of mobile component species simulated + icbndflg : int (default is 1) + a flag that determines if active domain for transport the same as for flow + itvd : int (default is 3) + 0 : upstream weighted scheme is used for simulating the advective term + >0 : number of TVD correction iterations used for simulating the advective term + iadsorb : int (default is 0) + adsorption flag (0: no adsorption, 1: linear isotherm, 2: Freundlich isotherm, + 3: Langmuir isotherm) + ict : int (default is 0) + transport solution scheme (0:water phase concentration, 1:total concentration) + cinact : float (default is -1.0) + concentration value that will be output at inactive nodes + ciclose : float (default is 1e-6) + concentration tolerance for convergence of the matrix solver + idisp : int (default is 1) + flag of dispersion (0: no dispersion, 1: isotropic, 2: anisotropic) + ixdisp : int (default is 0) + flag of cross-dispersion (0: no cross-dispersion, 1: cross-dispersion) + diffnc : float (default is 0.57024) + molecular diffusion coefficient + izod : int (default is 0) + flag of zero-order decay (0: no zero-order decay, 1: in water + 2: on soil, 3: on water and soil, 4: on air-water interface) + ifod : int (default is 0) + flag of first-order decay (0: no first-order decay, 1: in water + 2: on soil, 3: on water and soil, 4: on air-water interface) + ifmbc : int (default is 0) + flag of flux mass balance errors (0: not considered, 1: computed and reported) + iheat : int (default is 0) + flag of energy balance equation (0: not considered, 1: computed) + imcomp : int (default is 0) + number of immobile component species simulated + idispcln : int (default is 0) + index connection between GWF and CLN cells (0: finite difference approximation + 1: Thiem solution, 2: Thiem solution with local borehole thermal resistance) + nseqitr : int (default is 0) + an index or count for performing sequential iterations (0/1: no iterations + >1: number of iterations of the transport and reaction modules) + icbund : int or array of ints (nlay, nrow, ncol) + is the cell-by-cell flag for the transport simulation + prsity : float or array of floats (nlay, nrow, ncol), default is 0.15 + is the porosity of the porous medium + bulkd : float or array of floats (nlay, nrow, ncol), default is 157.0 + is the bulk density of the porous medium + anglex : float or array of floats (njag) + is the angle (in radians) between the horizontal x-axis and the outward + normal to the face between a node and its connecting nodes. The angle + varies between zero and 6.283185 (two pi being 360 degrees). + dl : float or array of floats (nlay, nrow, ncol), default is 1.0 + longitudinal dispersivity + dt : float or array of floats (nlay, nrow, ncol), default is 0.1 + transverse dispersivity + dlx : float or array of floats (nlay-1, nrow, ncol), default is 1.0 + x-direction longitudinal dispersivity + dly : float or array of floats (nlay, nrow, ncol), default is 1.0 + y-direction longitudinal dispersivity + dlz : float or array of floats (nlay, nrow, ncol), default is 0.1 + z-direction longitudinal dispersivity + dtxy : float or array of floats (nlay, nrow, ncol), default is 0.1 + xy-direction transverse dispersivity + dtyz : float or array of floats (nlay, nrow, ncol), default is 0.1 + yz-direction transverse dispersivity + dtxz : float or array of floats (nlay, nrow, ncol), default is 0.1 + xz-direction transverse dispersivity + adsorb : float or array of floats (nlay, nrow, ncol), default is none + adsorption coefficient of a contaminant species + flich : float or array of floats (nlay, nrow, ncol), default is none + Freundlich adsorption isotherm exponent + zodrw : float or array of floats (nlay, nrow, ncol), default is none + zero-order decay coefficient in water + zodrs : float or array of floats (nlay, nrow, ncol),default is none + zero-order decay coefficient on soil + zodraw : float or array of floats (nlay, nrow, ncol), default is none + zero-order decay coefficient on air-water interface + fodrw : float or array of floats (nlay, nrow, ncol), default is none + first-order decay coefficient in water + fodrs : float or array of floats (nlay, nrow, ncol), default is none + first-order decay coefficient on soil + fodraw : float or array of floats (nlay, nrow, ncol), default is none + first-order decay coefficient on air-water interface + conc : float or array of floats (nlay, nrow, ncol), default is 0.0 + initial concentration of each contaminant species + extension : string + mbegwunf : flow imbalance information + mbegwunt : transport mass imbalance information + mbeclnunf : flow imbalance information for CLN domain + mbeclnunt : transport mass imbalance information for the CLN domain + (default is ['bct','mbegwunf','mbegwunt','mbeclnunf','mbeclnunt']). + unitnumber : int + File unit number and the output files. + (default is [59, 0, 0, 0, 0, 0, 0] ). + filenames : str or list of str + Filenames to use for the package and the output files. + add_package : bool, default is True + Flag to add the initialised package object to the parent model object. + + Methods + ------- + + See Also + -------- + + Notes + ----- + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg(exe_name='USGs_1.exe') + >>> disu = flopy.mfusg.MfUsgDisU(model=ml, nlay=1, nodes=1, + iac=[1], njag=1,ja=np.array([0]), fahl=[1.0], cl12=[1.0]) + >>> bct = flopy.mfusg.MfUsgBct(ml)""" + + def __init__( + self, + model, + itrnsp=1, + ipakcb=0, + mcomp=1, + icbndflg=1, + itvd=3, + iadsorb=0, + ict=0, + cinact=-1.0, + ciclose=1e-6, + idisp=1, + ixdisp=0, + diffnc=0.57024, + izod=0, + ifod=0, + ifmbc=0, + iheat=0, + imcomp=0, + idispcln=0, + nseqitr=0, + timeweight=None, + chaindecay=False, + nparent=0, + jparent=0, + stotio=0, + sptlrct=0, + only_satadsorb=False, + spatialreact=False, + solubility=False, + aw_adsorb=False, + iarea_fn=None, + ikawi_fn=None, + imasswr=None, + crootname=None, + mbegwunf=0, + mbegwunt=0, + mbeclnunf=0, + mbeclnunt=0, + htcondw=0, + rhow=0, + htcapw=0, + htcaps=0, + htconds=0, + heat=0.0, + icbund=1, + prsity=0.15, + bulkd=157.0, + anglex=None, + dl=1.0, + dt=0.1, + dlx=1.0, + dly=1.0, + dlz=0.1, + dtxy=0.1, + dtyz=0.1, + dtxz=0.1, + sollim=0.0, + solslope=0.0, + alangaw=0.0, + blangaw=0.0, + adsorb=0, + flich=0, + zodrw=0, + zodrs=0, + zodraw=0, + fodrw=0, + fodrs=0, + fodraw=0, + conc=0.0, + imconc=0.0, + extension=["bct", "cbt", "mbegwf", "mbegwt", "mbeclnf", "mbeclnt"], + unitnumber=None, + filenames=None, + add_package=True, + ): + """Constructs the MfUsgBct object.""" + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if unitnumber is None: + unitnumber = MfUsgBct._defaultunit() + elif isinstance(unitnumber, list): + if len(unitnumber) < 6: + for idx in range(len(unitnumber), 6): + unitnumber.append(0) + + # set filenames + filenames = self._prepare_filenames(filenames, 6) + + # cbc output file + self.set_cbc_output_file(ipakcb, model, filenames[1]) + + super().__init__( + model, + extension, + self._ftype(), + unitnumber, + filenames, + ) + + self._generate_heading() + self.itrnsp = itrnsp + self.ipakcb = ipakcb + self.mcomp = mcomp + self.icbndflg = icbndflg + self.itvd = itvd + self.iadsorb = iadsorb + self.ict = ict + self.cinact = cinact + self.ciclose = ciclose + self.idisp = idisp + self.ixdisp = ixdisp + self.diffnc = diffnc + self.izod = izod + self.ifod = ifod + self.ifmbc = ifmbc + self.iheat = iheat + self.imcomp = imcomp + self.idispcln = idispcln + self.nseqitr = nseqitr + + model.mcomp = mcomp + model.iheat = iheat + + self.timeweight = timeweight + self.chaindecay = chaindecay + self.only_satadsorb = only_satadsorb + self.spatialreact = spatialreact + self.solubility = solubility + self.aw_adsorb = aw_adsorb + self.iarea_fn = iarea_fn + self.ikawi_fn = ikawi_fn + self.imasswr = imasswr + self.crootname = crootname + + self.mbegwunf = mbegwunf + self.mbegwunt = mbegwunt + self.mbeclnunf = mbeclnunf + self.mbeclnunt = mbeclnunt + + self.htcondw = htcondw + self.rhow = rhow + self.htcapw = htcapw + + nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper + + if self.icbndflg == 0: + if icbund is None: + icbund = 1 + self.icbund = Util3d( + model, (nlay, nrow, ncol), np.int32, icbund, name="icbund" + ) + + self.prsity = Util3d( + model, (nlay, nrow, ncol), np.float32, prsity, name="prsity" + ) + + if self.iadsorb or self.iheat: + if bulkd is None: + bulkd = 157.0 + self.bulkd = Util3d( + model, (nlay, nrow, ncol), np.float32, bulkd, name="bulkd" + ) + + if anglex is not None: + njag = self.parent.get_package("DISU").njag + self.anglex = Util2d(model, (njag,), np.float32, anglex, name="anglex") + + if self.idisp == 1: + self.dl = Util3d(model, (nlay, nrow, ncol), np.float32, dl, name="dl") + self.dt = Util3d(model, (nlay, nrow, ncol), np.float32, dt, name="dt") + + if self.idisp == 2: + self.dlx = Util3d(model, (nlay, nrow, ncol), np.float32, dlx, name="dlx") + self.dly = Util3d(model, (nlay, nrow, ncol), np.float32, dly, name="dly") + self.dlz = Util3d(model, (nlay, nrow, ncol), np.float32, dlz, name="dlz") + self.dtxy = Util3d(model, (nlay, nrow, ncol), np.float32, dtxy, name="dtxy") + self.dtyz = Util3d(model, (nlay, nrow, ncol), np.float32, dtyz, name="dtyz") + self.dtxz = Util3d(model, (nlay, nrow, ncol), np.float32, dtxz, name="dtxz") + if self.iheat == 1: + self.htcaps = Util3d( + model, (nlay, nrow, ncol), np.float32, htcaps, name="htcaps" + ) + self.htconds = Util3d( + model, (nlay, nrow, ncol), np.float32, htconds, name="htconds" + ) + if self.solubility == 1: + self.sollim = Util2d(model, (mcomp,), np.float32, sollim, name="sollim") + self.solslope = Util2d( + model, (mcomp,), np.float32, solslope, name="solslope" + ) + + if isinstance(adsorb, (int, float)): + adsorb = [adsorb] * mcomp + if isinstance(flich, (int, float)): + flich = [flich] * mcomp + if isinstance(zodrw, (int, float)): + zodrw = [zodrw] * mcomp + if isinstance(zodrs, (int, float)): + zodrs = [zodrs] * mcomp + if isinstance(zodraw, (int, float)): + zodraw = [zodraw] * mcomp + if isinstance(fodrw, (int, float)): + fodrw = [fodrw] * mcomp + if isinstance(fodrs, (int, float)): + fodrs = [fodrs] * mcomp + if isinstance(fodraw, (int, float)): + fodraw = [fodraw] * mcomp + if isinstance(conc, (int, float)): + conc = [conc] * mcomp + + self.adsorb = [0] * mcomp + self.flich = [0] * mcomp + self.zodrw = [0] * mcomp + self.zodrs = [0] * mcomp + self.zodraw = [0] * mcomp + self.fodrw = [0] * mcomp + self.fodrs = [0] * mcomp + self.fodraw = [0] * mcomp + self.conc = [0] * mcomp + + self.nparent = nparent + self.jparent = [0] * mcomp + self.stotio = [0] * mcomp + self.sptlrct = [0] * mcomp + + if isinstance(alangaw, (int, float)): + alangaw = [alangaw] * mcomp + if isinstance(blangaw, (int, float)): + blangaw = [blangaw] * mcomp + self.alangaw = [0] * mcomp + self.blangaw = [0] * mcomp + + for icomp in range(mcomp): + if self.chaindecay and self.nparent[icomp] > 0: + self.jparent[icomp] = Util2d( + model, (mcomp,), np.float32, jparent[icomp], name="jparent" + ) + self.stotio[icomp] = Util2d( + model, (mcomp,), np.float32, stotio[icomp], name="stotio" + ) + # Not tested + if self.spatialreact: + self.sptlrct[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + sptlrct[icomp], + name="sptlrct", + ) + + if self.aw_adsorb and self.iarea_fn in (1, 2, 3): + self.alangaw[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + alangaw[icomp], + name="alangaw species {icomp+1}", + ) + self.blangaw[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + blangaw[icomp], + name="blangaw species {icomp+1}", + ) + + self.adsorb[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + adsorb[icomp], + name=f"adsorb species {icomp+1}", + ) + self.flich[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + flich[icomp], + name=f"flich species {icomp+1}", + ) + self.zodrw[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + zodrw[icomp], + name=f"zodrw species {icomp+1}", + ) + self.zodrs[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + zodrs[icomp], + name=f"zodrs species {icomp+1}", + ) + self.zodraw[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + zodraw[icomp], + name=f"zodraw species {icomp+1}", + ) + self.fodrw[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + fodrw[icomp], + name=f"fodrw species {icomp+1}", + ) + self.fodrs[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + fodrs[icomp], + name=f"fodrs species {icomp+1}", + ) + self.fodraw[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + fodraw[icomp], + name=f"fodraw species {icomp+1}", + ) + self.conc[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + conc[icomp], + name=f"conc species {icomp+1}", + ) + + if self.iheat == 1: + self.heat = Util3d(model, (nlay, nrow, ncol), np.float32, heat, name="heat") + + if self.imcomp > 0: + if isinstance(imconc, (int, float)): + imconc = [imconc] * mcomp + self.imconc = [0] * mcomp + for icomp in range(self.imcomp): + self.imconc[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + imconc[icomp], + name=f"imconc species {icomp+1}", + ) + + if add_package: + self.parent.add_package(self) + + def write_file(self, f=None): + """ + Write the BCT package file. + + Parameters + ---------- + f : open file object. + Default is None, which will result in MfUsg.fn_path being + opened for writing. + + """ + # Open file for writing + if f is None: + f_obj = open(self.fn_path, "w") + + # Item 1: ITRNSP ipakcb MCOMP ICBNDFLG ITVD IADSORB ICT CINACT CICLOSE + # IDISP IXDISP DIFFNC IZOD IFOD IFMBC IHEAT IMCOMP IDISPCLN NSEQITR + f_obj.write(f"{self.heading}\n") + f_obj.write( + f" {self.itrnsp:9d} {self.ipakcb:9d} {self.mcomp:9d} {self.icbndflg:9d} " + f"{self.itvd:9d} {self.iadsorb:9d} {self.ict:9d} {self.cinact:9.2f} " + f"{self.ciclose:9.2e} {self.idisp:9d} {self.ixdisp:9d} {self.diffnc:9.2f} " + f"{self.izod:9d} {self.ifod:9d} {self.ifmbc:9d} {self.iheat:9d} " + f"{self.imcomp:9d} {self.idispcln:9d} {self.nseqitr:9d}" + ) + + # Options: TIMEWEIGHT CHAINDECAY ONLY_SATADSORB SPATIALREACT SOLUBILITY + # A-W_ADSORB WRITE_GWMASS MULTIFILE + if self.chaindecay: + f_obj.write(" CHAINDECAY") + if self.only_satadsorb: + f_obj.write(" ONLY_SATADSORB") + if self.spatialreact: + f_obj.write(" SPATIALREACT") + if self.solubility: + f_obj.write(" SOLUBILITY") + if self.timeweight is not None: + f_obj.write(f" TIMEWEIGHT {self.timeweight:9.2f}") + if self.aw_adsorb: + f_obj.write(f" A-W_ADSORB {self.iarea_fn:9d} {self.ikawi_fn:9d}") + if self.imasswr is not None: + f_obj.write(f" WRITE_GWMASS {self.imasswr:9d}") + if self.crootname is not None: + f_obj.write(f" MULTIFILE {self.crootname}") + + f_obj.write("\n") + + # Item 1b: IFMBC == 1 + if self.ifmbc == 1: + f_obj.write( + f" {self.unit_number[1]:9d} {self.unit_number[2]:9d}" + f" {self.unit_number[3]:9d} {self.unit_number[4]:9d}\n" + ) + + # Item 1c: IHEAT == 1 + if self.iheat == 1: + f_obj.write(f" {self.htcondw:9.2e} {self.rhow:9.2e} {self.htcapw:9.2e}\n") + + # Not implemented - item 1d aw_adsorb used and iarea_fn = 5 or ikawi_fn = 4 + # "nazones": int + # "natabrows": int + # Not implemented - item 1e aw_adsorb used and iarea_fn = 5 or ikawi_fn = 4 + # "iawizonmap": int + # Not implemented - item 1f aw_adsorb used and iarea_fn = 3 + # "rog_sigma": int + # Not implemented - item 1g aw_adsorb used and ikawi_fn = 3 + # "sigma_rt": int + # Not implemented - item 1h aw_adsorb used and iarea_fn = 5 + # "awi_area_tab": int + + # Item 2: ICBUND + if self.icbndflg == 0: + f_obj.write(self.icbund.get_file_entry()) + + # Item 3: PRSITY + f_obj.write(self.prsity.get_file_entry()) + + # Item 4: BULKD + if self.iadsorb or self.iheat: + f_obj.write(self.bulkd.get_file_entry()) + + # Item 5: ANGLEX + structured = self.parent.structured + if (not structured) and self.idisp: + f_obj.write(self.anglex.get_file_entry()) + + # Item 6 & 7: DL & DT + if self.idisp == 1: + f_obj.write(self.dl.get_file_entry()) + f_obj.write(self.dt.get_file_entry()) + + # Item 8 - 13: DLX, DLY, DLZ, DTXY, DTYZ, DTXZ + if self.idisp == 2: + f_obj.write(self.dlx.get_file_entry()) + f_obj.write(self.dly.get_file_entry()) + f_obj.write(self.dlz.get_file_entry()) + f_obj.write(self.dtxy.get_file_entry()) + f_obj.write(self.dtyz.get_file_entry()) + f_obj.write(self.dtxz.get_file_entry()) + + if self.iheat == 1: + f_obj.write(self.htcaps.get_file_entry()) + f_obj.write(self.htconds.get_file_entry()) + + for icomp in range(self.mcomp): + if self.chaindecay: + f_obj.write(f"{self.nparent[icomp]:9d}\n") + if self.nparent[icomp] > 0: + f_obj.write(self.jparent[icomp].get_file_entry()) + f_obj.write(self.stotio[icomp].get_file_entry()) + if self.spatialreact: + f_obj.write(self.sptlrct[icomp].get_file_entry()) + + # Item A6-A7: ALANGAW, BLANGAW + if self.aw_adsorb and self.iarea_fn in (1, 2, 3): + f_obj.write(self.alangaw[icomp].get_file_entry()) + f_obj.write(self.blangaw[icomp].get_file_entry()) + + # Item 14 - 19: ADSORB, FLICH, ZODRW, ZODRS, ZODRAW, FODRW, FODRS, FODRAW + if self.iadsorb: + f_obj.write(self.adsorb[icomp].get_file_entry()) + if self.iadsorb == 2 or self.iadsorb == 3: + f_obj.write(self.flich[icomp].get_file_entry()) + if self.izod == 1 or self.izod == 3 or self.izod == 4: + f_obj.write(self.zodrw[icomp].get_file_entry()) + if self.iadsorb and (self.izod == 2 or self.izod == 3 or self.izod == 4): + f_obj.write(self.zodrs[icomp].get_file_entry()) + if self.aw_adsorb and self.izod == 4: + f_obj.write(self.zodraw[icomp].get_file_entry()) + if self.ifod == 1 or self.ifod == 3 or self.ifod == 4: + f_obj.write(self.fodrw[icomp].get_file_entry()) + if self.iadsorb and (self.ifod == 2 or self.ifod == 3 or self.ifod == 4): + f_obj.write(self.fodrs[icomp].get_file_entry()) + if self.aw_adsorb and self.ifod == 4: + f_obj.write(self.fodraw[icomp].get_file_entry()) + + # Item 20: CONC + f_obj.write(self.conc[icomp].get_file_entry()) + + if self.iheat == 1: + f_obj.write(self.heat.get_file_entry()) + + if self.imcomp > 0: + for icomp in range(self.imcomp): + f_obj.write(self.imconc[icomp].get_file_entry()) + + # close the file + f_obj.close() + + @classmethod + def load(cls, f, model, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + bct : MfUsgBct object + + Examples + -------- + + Examples + -------- + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> dis = flopy.modflow.ModflowDis.load('Test1.dis', ml) + >>> bct = flopy.mfusg.MfUsgBct.load('Test1.BTN', ml) + + """ + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + # determine problem dimensions + nlay = model.nlay + dis = model.get_package("DIS") + if dis is None: + dis = model.get_package("DISU") + njag = dis.njag + + if model.verbose: + print("loading bct package file...") + + f_obj = get_open_file_object(f, "r") + + # item 0 + line = f_obj.readline().upper() + while line.startswith("#"): + line = f_obj.readline().upper() + + t = line.split() + kwargs = {} + + # item 1a + vars = { + "itrnsp": int, + "ipakcb": int, + "mcomp": int, + "icbndflg": int, + "itvd": int, + "iadsorb": int, + "ict": int, + "cinact": float, + "ciclose": float, + "idisp": int, + "ixdisp": int, + "diffnc": float, + } + + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = c(t[i].strip()) + print(f"{v}={kwargs[v]}") + + vars = { + "izod": int, + "ifod": int, + "ifmbc": int, + "iheat": int, + "imcomp": int, + "idispcln": int, + "nseqitr": int, + } + + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = type_from_iterable(t, index=12 + i, _type=c, default_val=0) + print(f"{v}={kwargs[v]}") + + # item 1a - options + if "TIMEWEIGHT" in t: + idx = t.index("TIMEWEIGHT") + kwargs["timeweight"] = float(t[idx + 1]) + else: + kwargs["timeweight"] = None + + if "CHAINDECAY" in t: + kwargs["chaindecay"] = 1 + else: + kwargs["chaindecay"] = 0 + + if "ONLY_SATADSORB" in t: + kwargs["only_satadsorb"] = 1 + else: + kwargs["only_satadsorb"] = 0 + + if "SPATIALREACT" in t: + kwargs["spatialreact"] = 1 + else: + kwargs["spatialreact"] = 0 + + if "SOLUBILITY" in t: + kwargs["solubility"] = 1 + else: + kwargs["solubility"] = 0 + + if "A-W_ADSORB" in t: + idx = t.index("A-W_ADSORB") + kwargs["aw_adsorb"] = 1 + kwargs["iarea_fn"] = int(t[idx + 1]) + kwargs["ikawi_fn"] = int(t[idx + 2]) + else: + kwargs["aw_adsorb"] = 0 + kwargs["iarea_fn"] = None + kwargs["ikawi_fn"] = None + + if "WRITE_GWMASS" in t: + idx = t.index("WRITE_GWMASS") + kwargs["imasswr"] = int(t[idx + 1]) + else: + kwargs["imasswr"] = None + + if "MULTIFILE" in t: + idx = t.index("MULTIFILE") + kwargs["crootname"] = str(t[idx + 1]) + else: + kwargs["crootname"] = None + + # item 1b ifmbc == 1 + if kwargs["ifmbc"] == 1: + vars = { + "mbegwunf": int, + "mbegwunt": int, + "mbeclnunf": int, + "mbeclnunt": int, + } + line = f_obj.readline().upper() + t = line.split() + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = c(t[i].strip()) + else: + kwargs["mbegwunf"] = 0 + kwargs["mbegwunt"] = 0 + kwargs["mbeclnunf"] = 0 + kwargs["mbeclnunt"] = 0 + + # item 1c iheat == 1 + if kwargs["iheat"]: + vars = { + "htcondw": float, + "rhow": float, + "htcapw": float, + } + line = f_obj.readline().upper() + t = line.split() + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = c(t[i].strip()) + + # Not implemented - item 1d aw_adsorb used and iarea_fn = 5 or ikawi_fn = 4 + # "nazones": int + # "natabrows": int + # Not implemented - item 1e aw_adsorb used and iarea_fn = 5 or ikawi_fn = 4 + # "iawizonmap": int + # Not implemented - item 1f aw_adsorb used and iarea_fn = 3 + # "rog_sigma": int + # Not implemented - item 1g aw_adsorb used and ikawi_fn = 3 + # "sigma_rt": int + # Not implemented - item 1h aw_adsorb used and iarea_fn = 5 + # "awi_area_tab": int + + # item 2 ICBUND + if kwargs["icbndflg"] == 0: + kwargs["icbund"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "icbund", ext_unit_dict + ) + + # item 3 PRSITY + kwargs["prsity"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "prsity", ext_unit_dict + ) + + # item 4 BULKD + if kwargs["iadsorb"] or kwargs["iheat"]: + kwargs["bulkd"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "bulkd", ext_unit_dict + ) + + # item 5 ANGLEX + if (not model.structured) and kwargs["idisp"]: + if model.verbose: + print(" loading ANGLEX...") + kwargs["anglex"] = Util2d.load( + f_obj, model, (njag,), np.float32, "anglex", ext_unit_dict + ) + + # item 6 & 7 DL & DT + if kwargs["idisp"] == 1: + kwargs["dl"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dl", ext_unit_dict + ) + kwargs["dt"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dt", ext_unit_dict + ) + + # item 8 - 13 DLX, DLY, DLZ, DTXY, DTYZ, DTXZ + if kwargs["idisp"] == 2: + kwargs["dlx"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dlx", ext_unit_dict + ) + kwargs["dly"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dly", ext_unit_dict + ) + kwargs["dlz"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dlz", ext_unit_dict + ) + kwargs["dtxy"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dtxy", ext_unit_dict + ) + kwargs["dtyz"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dtyz", ext_unit_dict + ) + kwargs["dtxz"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dtxz", ext_unit_dict + ) + + # item H1 AND H2 + if kwargs["iheat"] == 1: + kwargs["htcaps"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "htcondw", ext_unit_dict + ) + kwargs["htconds"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "rhow", ext_unit_dict + ) + + # Not implemented item A1 - A5 if AW_ADSORB option is on + + mcomp = kwargs["mcomp"] + + adsorb = [0] * mcomp + flich = [0] * mcomp + zodrw = [0] * mcomp + zodrs = [0] * mcomp + zodraw = [0] * mcomp + fodrw = [0] * mcomp + fodrs = [0] * mcomp + fodraw = [0] * mcomp + conc = [0] * mcomp + + nparent = np.zeros(mcomp, dtype=np.int32) + jparent = [0] * mcomp + stotio = [0] * mcomp + sptlrct = [0] * mcomp + + sollim = [0] * mcomp + solslope = [0] * mcomp + + alangaw = [0] * mcomp + blangaw = [0] * mcomp + + for icomp in range(mcomp): + # item S1 - S2 if SOLUBILITY option is on -- not tested + if kwargs["solubility"] == 1: + if model.verbose: + print(" loading SOLUBILITY ...") + sollim[icomp] = read1d(f_obj, sollim) + solslope[icomp] = read1d(f_obj, solslope) + + # item C1 - C3 if CHAINDECAY option is on + if kwargs["chaindecay"] == 1: + if model.verbose: + print(" loading CHAINDECAY...") + t = f_obj.readline().split() + # print(t) + nparent[icomp] = int(t[0].strip()) + if nparent[icomp] > 0: + jparent[icomp] = Util2d.load( + f_obj, model, (mcomp,), np.int32, "jparent", ext_unit_dict + ) + stotio[icomp] = Util2d.load( + f_obj, model, (mcomp,), np.float32, "stotio", ext_unit_dict + ) + # item C4 if SPATIALREACT option is on - Not tested + if kwargs["spatialreact"] == 1: + if model.verbose: + print(" loading SPATIALREACT...") + sptlrct[icomp] = Util3d.load( + f_obj, model, (nlay,), np.float32, "sptlrct", ext_unit_dict + ) + + # item A6 - A7 if AW_ADSORB option is on + if kwargs["aw_adsorb"] == 1 and kwargs["iarea_fn"] in (1, 2, 3): + alangaw[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "alangaw", ext_unit_dict + ) + blangaw[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "blangaw", ext_unit_dict + ) + + # item A8 not implemented if AW_ADSORB option is on and iarea_fn = 4 + + # item 14 - 20 ADSORB, FLICH, ZODRW, ZODRS, ZODRAW, FODRW, FODRS, FODRAW + if kwargs["iadsorb"]: + adsorb[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "adsorb", ext_unit_dict + ) + if kwargs["iadsorb"] == 2 or kwargs["iadsorb"] == 3: + flich[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "flich", ext_unit_dict + ) + if kwargs["izod"] == 1 or kwargs["izod"] == 3 or kwargs["izod"] == 4: + zodrw[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "zodrw", ext_unit_dict + ) + if kwargs["iadsorb"] and ( + kwargs["izod"] == 2 or kwargs["izod"] == 3 or kwargs["izod"] == 4 + ): + zodrs[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "zodrs", ext_unit_dict + ) + if kwargs["aw_adsorb"] and (kwargs["izod"] == 4): + zodraw[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "zodraw", ext_unit_dict + ) + if kwargs["ifod"] == 1 or kwargs["ifod"] == 3 or kwargs["ifod"] == 4: + fodrw[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "fodrw", ext_unit_dict + ) + if kwargs["iadsorb"] and ( + kwargs["ifod"] == 2 or kwargs["ifod"] == 3 or kwargs["ifod"] == 4 + ): + fodrs[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "fodrs", ext_unit_dict + ) + if kwargs["aw_adsorb"] and (kwargs["ifod"] == 4): + fodraw[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "fodraw", ext_unit_dict + ) + + # item 20 CONC + conc[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "conc", ext_unit_dict + ) + + if kwargs["iheat"] == 1: + kwargs["heat"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "heat", ext_unit_dict + ) + + if kwargs["imcomp"] > 0: + imconc = [0] * kwargs["imcomp"] + for icomp in range(kwargs["imcomp"]): + imconc[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "imconc", ext_unit_dict + ) + kwargs["imconc"] = imconc + + kwargs["nparent"] = nparent + kwargs["jparent"] = jparent + kwargs["stotio"] = stotio + kwargs["sptlrct"] = sptlrct + kwargs["sollim"] = sollim + kwargs["solslope"] = solslope + kwargs["alangaw"] = alangaw + kwargs["blangaw"] = blangaw + + kwargs["adsorb"] = adsorb + kwargs["flich"] = flich + kwargs["zodrw"] = zodrw + kwargs["zodrs"] = zodrs + kwargs["zodraw"] = zodraw + kwargs["fodrw"] = fodrw + kwargs["fodrs"] = fodrs + kwargs["fodraw"] = fodraw + kwargs["conc"] = conc + + f_obj.close() + + # set package unit number + # reset unit numbers + unitnumber = cls._defaultunit() + filenames = [None] * 6 + unitnumber[0], filenames[0] = model.get_ext_dict_attr( + ext_unit_dict, filetype=cls._ftype() + ) + if kwargs["ipakcb"] > 0: + unitnumber[1], filenames[1] = model.get_ext_dict_attr( + ext_unit_dict, unit=kwargs["ipakcb"] + ) + if kwargs["mbegwunf"] > 0: + unitnumber[2], filenames[2] = model.get_ext_dict_attr( + ext_unit_dict, unit=kwargs["mbegwunf"] + ) + if kwargs["mbegwunt"] > 0: + unitnumber[3], filenames[3] = model.get_ext_dict_attr( + ext_unit_dict, unit=kwargs["mbegwunt"] + ) + if kwargs["mbeclnunf"] > 0: + unitnumber[4], filenames[4] = model.get_ext_dict_attr( + ext_unit_dict, unit=kwargs["mbeclnunf"] + ) + if kwargs["mbeclnunt"] > 0: + unitnumber[5], filenames[5] = model.get_ext_dict_attr( + ext_unit_dict, unit=kwargs["mbeclnunt"] + ) + + return cls(model, unitnumber=unitnumber, filenames=filenames, **kwargs) + + @staticmethod + def _load_prop_arrays(f_obj, model, nlay, dtype, name, ext_unit_dict): + if model.verbose: + print(f" loading {name} ...") + prop_array = [0] * nlay + for layer in range(nlay): + util2d_shape = get_util2d_shape_for_layer(model, layer=layer) + prop_array[layer] = Util2d.load( + f_obj, model, util2d_shape, dtype, name, ext_unit_dict + ) + return prop_array + + @staticmethod + def _ftype(): + return "BCT" + + @staticmethod + def _defaultunit(): + return [150, 0, 0, 0, 0, 0] diff --git a/flopy/mfusg/mfusgcln.py b/flopy/mfusg/mfusgcln.py index e88b8f312c..d5ba53015e 100644 --- a/flopy/mfusg/mfusgcln.py +++ b/flopy/mfusg/mfusgcln.py @@ -196,6 +196,7 @@ def __init__( ), unitnumber=None, filenames=None, + **kwargs, ): """Package constructor.""" msg = ( @@ -311,6 +312,98 @@ def __init__( locat=self.unit_number[0], ) + # Transport parameters + bct = model.get_package("BCT") + if bct is not None: + if bct.icbndflg == 0: + icbund = kwargs.pop("icbund", None) + self.icbund = Util2d( + model, + (self.nclnnds,), + np.int32, + icbund, + name="icbund", + locat=self.unit_number[0], + ) + if bct.idisp: + dll = kwargs.pop("dll", None) + self.dll = Util2d( + model, + (self.nclnnds,), + np.float32, + dll, + name="dll", + locat=self.unit_number[0], + ) + dlm = kwargs.pop("dlm", None) + self.dlm = Util2d( + model, + (self.nclnnds,), + np.float32, + dlm, + name="dlm", + locat=self.unit_number[0], + ) + + mcomp = bct.mcomp + self.sptlrct = [0] * mcomp + self.zodrw = [0] * mcomp + self.fodrw = [0] * mcomp + self.conc = [0] * mcomp + + if bct.spatialreact: + sptlrct = kwargs.pop("sptlrct", None) + if isinstance(sptlrct, (int, float)): + sptlrct = [sptlrct] * mcomp + if bct.izod: + zodrw = kwargs.pop("zodrw", None) + if isinstance(zodrw, (int, float)): + zodrw = [zodrw] * mcomp + if bct.ifod: + fodrw = kwargs.pop("fodrw", None) + if isinstance(fodrw, (int, float)): + fodrw = [fodrw] * mcomp + conc = kwargs.pop("conc", None) + if isinstance(conc, (int, float)): + conc = [conc] * mcomp + + for icomp in range(mcomp): + if bct.spatialreact: + self.sptlrct[icomp] = Util2d( + model, + (self.nclnnds,), + np.float32, + sptlrct[icomp], + name="sptlrct", + locat=self.unit_number[0], + ) + if bct.izod: + self.zodrw[icomp] = Util2d( + model, + (self.nclnnds,), + np.float32, + zodrw[icomp], + name="zodrw", + locat=self.unit_number[0], + ) + if bct.ifod: + self.fodrw[icomp] = Util2d( + model, + (self.nclnnds,), + np.float32, + fodrw[icomp], + name="fodrw", + locat=self.unit_number[0], + ) + self.conc[icomp] = Util2d( + model, + (self.nclnnds,), + np.float32, + conc[icomp], + name="conc", + locat=self.unit_number[0], + ) + self.parent.add_package(self) @staticmethod @@ -484,6 +577,8 @@ def write_file(self, f=None, check=False): f_cln.write(self.ibound.get_file_entry()) f_cln.write(self.strt.get_file_entry()) + self._write_transport(f_cln) + f_cln.close() def _write_items_0_1(self, f_cln): @@ -503,19 +598,37 @@ def _write_items_0_1(self, f_cln): ) if self.nrectyp > 0: - f_cln.write(f"RECTANGULAR {self.nrectyp:d}") + f_cln.write(f" RECTANGULAR {self.nrectyp:d}") if self.bhe: - f_cln.write("BHEDETAIL ") + f_cln.write(" BHEDETAIL ") if self.iclncn != 0: - f_cln.write(f"SAVECLNCON {self.iclncn:d}") + f_cln.write(f" SAVECLNCON {self.iclncn:d}") if self.iclnmb != 0: - f_cln.write(f"SAVECLNMAS {self.iclnmb:d}") + f_cln.write(f" SAVECLNMAS {self.iclnmb:d}") if self.grav is not None: - f_cln.write(f"GRAVITY {self.grav:f}") + f_cln.write(f" GRAVITY {self.grav:f}") if self.visk is not None: - f_cln.write(f"VISCOSITY {self.visk:f}") + f_cln.write(f" VISCOSITY {self.visk:f}") f_cln.write("\n") + def _write_transport(self, f_cln): + bct = self.parent.get_package("BCT") + + if bct is not None: + if bct.icbndflg == 0: + f_cln.write(self.icbund.get_file_entry()) + if bct.idisp: + f_cln.write(self.dll.get_file_entry()) + f_cln.write(self.dlm.get_file_entry()) + for icomp in range(bct.mcomp): + if bct.spatialreact: + f_cln.write(self.sptlrct[icomp].get_file_entry()) + if bct.izod: + f_cln.write(self.zodrw[icomp].get_file_entry()) + if bct.ifod: + f_cln.write(self.fodrw[icomp].get_file_entry()) + f_cln.write(self.conc[icomp].get_file_entry()) + @classmethod def load(cls, f, model, pak_type="cln", ext_unit_dict=None, **kwargs): """ @@ -615,6 +728,20 @@ def load(cls, f, model, pak_type="cln", ext_unit_dict=None, **kwargs): print(" Reading strt...") strt = Util2d.load(f, model, (nclnnds, 1), np.float32, "strt", ext_unit_dict) + bct = model.get_package("BCT") + if bct is not None: + if model.verbose: + print("loading transport parameters...") + ( + kwargs["icbund"], + kwargs["dll"], + kwargs["dlm"], + kwargs["sptlrct"], + kwargs["zodrw"], + kwargs["fodrw"], + kwargs["conc"], + ) = cls._load_transport(f, model, nclnnds, bct, ext_unit_dict) + if hasattr(f, "read"): f.close() @@ -662,6 +789,7 @@ def load(cls, f, model, pak_type="cln", ext_unit_dict=None, **kwargs): bhe=bhe, unitnumber=unitnumber, filenames=filenames, + **kwargs, ) # return dis object instance @@ -808,6 +936,51 @@ def _load_items_3to6(f_obj, model, ncln, iclnnds, ext_unit_dict): return nndcln, clncon, nja_cln, iac_cln, ja_cln, nclnnds + def _load_transport(f_obj, model, nclnnds, bct, ext_unit_dict): + """Loads cln items 3, or 4,5,6 from filehandle f.""" + + icbund = None + if bct.icbndflg == 0: + icbund = Util2d.load( + f_obj, model, (nclnnds, 1), np.int32, "icbund", ext_unit_dict + ) + + dll = None + dlm = None + if bct.idisp: + dll = Util2d.load( + f_obj, model, (nclnnds, 1), np.float32, "dll", ext_unit_dict + ) + dlm = Util2d.load( + f_obj, model, (nclnnds, 1), np.float32, "dlm", ext_unit_dict + ) + + mcomp = bct.mcomp + + sptlrct = [0] * mcomp + zodrw = [0] * mcomp + fodrw = [0] * mcomp + conc = [0] * mcomp + + for icomp in range(mcomp): + if bct.spatialreact: + sptlrct[icomp] = Util2d.load( + f_obj, model, (nclnnds, 1), np.float32, "sptlrct", ext_unit_dict + ) + if bct.izod: + zodrw[icomp] = Util2d.load( + f_obj, model, (nclnnds, 1), np.float32, "zodrw", ext_unit_dict + ) + if bct.ifod: + fodrw[icomp] = Util2d.load( + f_obj, model, (nclnnds, 1), np.float32, "fodrw", ext_unit_dict + ) + conc[icomp] = Util2d.load( + f_obj, model, (nclnnds, 1), np.float32, "conc", ext_unit_dict + ) + + return icbund, dll, dlm, sptlrct, zodrw, fodrw, conc + @staticmethod def _ftype(): return "CLN" diff --git a/flopy/mfusg/mfusgddf.py b/flopy/mfusg/mfusgddf.py new file mode 100644 index 0000000000..7cfd9ce8a9 --- /dev/null +++ b/flopy/mfusg/mfusgddf.py @@ -0,0 +1,219 @@ +""" +mfusgddf module. + +Contains the MfUsgDdf class. Note that the user can access +the MfUsgDdf class as `flopy.mfusg.MfUsgDdf`. +""" + +from ..pakbase import Package +from ..utils.flopy_io import line_parse +from ..utils.utils_def import ( + get_unitnumber_from_ext_unit_dict, + type_from_iterable, +) +from .mfusg import MfUsg + + +class MfUsgDdf(Package): + """MODFLOW-USG Transport Density Driven Flow (DDF) package class. + + parameters + ---------- + model : model object + the model object (of type :class:`flopy.mfusg.MfUsg`) to which + this package will be added. + rhofresh : density of freshwater + rhostd : density of standard solution + cstd : concentration of standard solution + ithickav : a flag indicating if thickness weighted averaging for density term + 0 = arithmetic averaging; 1 = thickness weighting averaging + imphdd : a flag of hydraulic head term in the density formulation + 0 = explicit (on the right-hand side vector) symmetry of the matrix + 1 = implicit (on the left-hand side matrix) asymmetric matrix + extension : str, optional + file extension (default is 'ddf'). + unitnumber : int, optional + fortran unit number for this package (default is none). + filenames : str or list of str + filenames to use for the package. if filenames=none the package name + will be created using the model name and package extension. if a + single string is passed the package will be set to the string. + default is none. + + attributes + ---------- + + methods + ------- + + see also + -------- + + notes + ----- + + examples + -------- + + >>> import flopy + >>> m = flopy.mfusg.mfusg() + >>> ddf = flopy.mfusg.mfusgddf(m) + """ + + def __init__( + self, + model, + rhofresh=1000.0, + rhostd=1025.0, + cstd=35.0, + ithickav=1, + imphdd=0, + extension="ddf", + unitnumber=None, + filenames=None, + ): + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + # set default unit number of one is not specified + if unitnumber is None: + unitnumber = self._defaultunit() + + # call base package constructor + super().__init__( + model, + extension=extension, + name=self._ftype(), + unit_number=unitnumber, + filenames=self._prepare_filenames(filenames), + ) + + self._generate_heading() + self.rhofresh = rhofresh + self.rhostd = rhostd + self.cstd = cstd + self.ithickav = ithickav + self.imphdd = imphdd + + self.parent.add_package(self) + return + + def write_file(self): + """ + Write the package file. + + Returns + ------- + None + + """ + f = open(self.fn_path, "w") + f.write(f"{self.heading}\n") + f.write( + f" {self.rhofresh:9.2f} {self.rhostd:9.2f} {self.cstd:9.2f}" + f" {self.ithickav:9d} {self.imphdd:9d}\n" + ) + f.close() + + @classmethod + def load(cls, f, model, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + ddf : MfUsgDdf object + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> ddf = flopy.mfusg.MfUsgDdf.load('Henry.ddf', ml) + """ + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if model.verbose: + print("loading ddf package file...") + + if model.version != "mfusg": + print( + "Warning: model version was reset from '{}' to 'mfusg' " + "in order to load a DDF file".format(model.version) + ) + model.version = "mfusg" + + openfile = not hasattr(f, "read") + if openfile: + filename = f + f = open(filename, "r") + + # dataset 0 -- header + line = f.readline().upper() + while line.startswith("#"): + line = f.readline().upper() + + if model.verbose: + print(" loading RHOFRESH RHOSTD CSTD ITHICKAV IMPHDD...") + + ll = line_parse(line) + rhofresh = float(ll.pop(0)) + rhostd = float(ll.pop(0)) + cstd = float(ll.pop(0)) + ithickav = type_from_iterable(ll, index=3, _type=int, default_val=1) + imphdd = type_from_iterable(ll, index=4, _type=int, default_val=0) + if model.verbose: + print( + f" RHOFRESH {rhofresh} \n RHOSTD {rhostd} \n CSTD {cstd} \n" + f" ITHICKAV {ithickav} \n IMPHDD {imphdd}" + ) + + if openfile: + f.close() + + # set package unit number + unitnumber, filenames = get_unitnumber_from_ext_unit_dict( + model, + cls, + ext_unit_dict, + ) + + return cls( + model, + rhofresh=rhofresh, + rhostd=rhostd, + cstd=cstd, + ithickav=ithickav, + imphdd=imphdd, + unitnumber=unitnumber, + filenames=filenames, + ) + + @staticmethod + def _ftype(): + return "DDF" + + @staticmethod + def _defaultunit(): + return 152 diff --git a/flopy/mfusg/mfusgdisu.py b/flopy/mfusg/mfusgdisu.py index 0c284387d2..709073d2c7 100644 --- a/flopy/mfusg/mfusgdisu.py +++ b/flopy/mfusg/mfusgdisu.py @@ -251,6 +251,8 @@ def __init__( # Set values of all parameters self._generate_heading() + model.structured = False + self.nodes = nodes self.nlay = nlay self.njag = njag diff --git a/flopy/mfusg/mfusgdpf.py b/flopy/mfusg/mfusgdpf.py new file mode 100644 index 0000000000..fb7223a8ae --- /dev/null +++ b/flopy/mfusg/mfusgdpf.py @@ -0,0 +1,327 @@ +""" +Mfusgdpf module. + +Contains the MfUsgDpf class. Note that the user can +access the MfUsgDpf class as `flopy.mfusg.MfUsgDpf`. +""" + +import numpy as np + +from ..pakbase import Package +from ..utils import Util2d, Util3d +from ..utils.utils_def import ( + get_open_file_object, + get_unitnumber_from_ext_unit_dict, + get_util2d_shape_for_layer, +) +from .mfusg import MfUsg + + +class MfUsgDpf(Package): + """Dual Porosity Flow (dpf) Package Class for MODFLOW-USG Transport. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.modflow.Modflow`) to which + this package will be added. + ipakcb : int (0,1,-1), (default is 0) + a flag and a unit number >0 for cell-by-cell mass flux terms. + idpfhd : int, (default is 0) + a flag and a unit number >0 for immobile domain heads. + idpfdd : int, (default is 0) + a flag and a unit number >0 for immobile domain drawdown. + iuzontabim : int or array of ints, (default is 0) + soil type index for all groundwater flow nodes in the immobile domain + iboundim : int, (default is 0) + boundary variable for the immobile domain (<0 constant head, =0 no flow) + hnewim : float or array of floats + initial (starting) head in the immobile domain + phif : float or array of floats + porosity of the mobile domain. + ddftr : float or array of floats + dual domain flow transfer rate + sc1im : float or array of floats + specific storage of the immobile domain + sc2im : float or array of floats + specific yield or porosity of the immobile domain + alphaim : float or array of floats + van Genuchten alpha coefficient of the immobile domain + betaim : float or array of floats + van Genuchten beta coefficient of the immobile domain + srim : float or array of floats + van Genuchten sr coefficient of the immobile domain + brookim : float or array of floats + Brooks-Corey exponent for the relative permeability of the immobile domain + bpim : float or array of floats + Bubble point or air entry pressure head of the immobile domain + extension : string, (default is 'dpf'). + unitnumber : int, default is 57. + File unit number. + filenames : str or list of str + Filenames to use for the package and the output files. + add_package : bool, default is True + Flag to add the initialised package object to the parent model object. + + Methods + ------- + + See Also + -------- + + Notes + ----- + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> disu = flopy.mfusg.MfUsgDisU(model=ml, nlay=1, nodes=1, + iac=[1], njag=1,ja=np.array([0]), fahl=[1.0], cl12=[1.0]) + >>> dpf = flopy.mfusg.MfUsgdpf(ml)""" + + def __init__( + self, + model, + ipakcb=0, + idpfhd=0, + idpfdd=0, + # iuzontabim=0, + iboundim=0, + hnewim=0.0, + phif=0.0, + ddftr=0.0, + sc1im=0.0, + sc2im=0.0, + # alphaim=0.0, + # betaim=0.0, + # srim=0.0, + # brookim=0.0, + # bpim=0.0, + extension="dpf", + unitnumber=None, + filenames=None, + add_package=True, + ): + """Constructs the MfUsgdpf object.""" + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + # set default unit number of one is not specified + if unitnumber is None: + self.unitnumber = self._defaultunit() + + super().__init__( + model, + extension, + self._ftype(), + unitnumber, + self._prepare_filenames(filenames), + ) + + self._generate_heading() + self.ipakcb = ipakcb + self.idpfhd = idpfhd + self.idpfdd = idpfdd + + model.idpf = 0 + + nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper + + # if self.parent.tabrich: + # self.iuzontabim = Util3d( + # model, (nlay, nrow, ncol), np.int32, iuzontabim, name="iuzontabim" + # ) + + self.iboundim = Util3d( + model, (nlay, nrow, ncol), np.float32, iboundim, name="iboundim" + ) + + self.hnewim = Util3d( + model, (nlay, nrow, ncol), np.float32, hnewim, name="hnewim" + ) + + self.phif = Util3d(model, (nlay, nrow, ncol), np.float32, phif, name="phif") + + self.ddftr = Util3d(model, (nlay, nrow, ncol), np.float32, ddftr, name="ddftr") + + self.sc1im = Util3d(model, (nlay, nrow, ncol), np.float32, sc1im, name="sc1im") + + self.sc2im = Util3d(model, (nlay, nrow, ncol), np.float32, sc2im, name="sc2im") + + # self.alphaim = Util3d( + # model, (nlay, nrow, ncol), np.float32, alphaim, name="alphaim" + # ) + + # self.betaim = Util3d( + # model, (nlay, nrow, ncol), np.float32, betaim, name="betaim" + # ) + + # self.srim = Util3d( + # model, (nlay, nrow, ncol), np.float32, srim, name="srim" + # ) + + # self.brookim = Util3d( + # model, (nlay, nrow, ncol), np.float32, brookim, name="brookim" + # ) + + # self.bpim = Util3d( + # model, (nlay, nrow, ncol), np.float32, bpim, name="bpim" + # ) + + if add_package: + self.parent.add_package(self) + + def write_file(self, f=None): + """ + Write the dpf package file. + + Parameters + ---------- + f : open file object. + Default is None, which will result in MfUsg.fn_path being + opened for writing. + + """ + # Open file for writing + if f is None: + f_obj = open(self.fn_path, "w") + + # f_obj.write(f"{self.heading}\n") + + # Item 0: IPAKCB, IdpfCON + f_obj.write(f" {self.ipakcb:9d} {self.idpfhd:9d} {self.idpfdd:9d} \n") + + f_obj.write(self.iboundim.get_file_entry()) + f_obj.write(self.hnewim.get_file_entry()) + f_obj.write(self.phif.get_file_entry()) + f_obj.write(self.ddftr.get_file_entry()) + f_obj.write(self.sc1im.get_file_entry()) + f_obj.write(self.sc2im.get_file_entry()) + + # close the file + f_obj.close() + + @classmethod + def load(cls, f, model, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + dpf : MfUsgdpf object + + Examples + -------- + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> dis = flopy.modflow.ModflowDis.load('Test1.dis', ml) + >>> dpf = flopy.mfusg.MfUsgdpf.load('Test1.BTN', ml) + """ + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if model.verbose: + print("loading dpf package file...") + + nlay = model.nlay + + f_obj = get_open_file_object(f, "r") + + # item 0 + line = f_obj.readline().upper() + while line.startswith("#"): + line = f_obj.readline().upper() + + t = line.split() + kwargs = {} + + # item 1a + vars = { + "ipakcb": int, + "idpfhd": int, + "idpfdd": int, + } + + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = c(t[i].strip()) + # print(f"{v}={kwargs[v]}\n") + + # item 1b + # if self.parent.tabrich: + # kwargs["iuzontabim"] = cls._load_prop_arrays( + # f_obj, model, nlay, np.int32, "iuzontabim", ext_unit_dict + # ) + + kwargs["iboundim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.int32, "iboundim", ext_unit_dict + ) + + kwargs["hnewim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "hnewim", ext_unit_dict + ) + + kwargs["phif"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "phif", ext_unit_dict + ) + + kwargs["ddftr"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "ddftr", ext_unit_dict + ) + + kwargs["sc1im"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "sc1im", ext_unit_dict + ) + + kwargs["sc2im"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "sc2im", ext_unit_dict + ) + + f_obj.close() + # set package unit number + unitnumber, filenames = get_unitnumber_from_ext_unit_dict( + model, cls, ext_unit_dict, kwargs["ipakcb"] + ) + + return cls(model, unitnumber=unitnumber, filenames=filenames, **kwargs) + + @staticmethod + def _load_prop_arrays(f_obj, model, nlay, dtype, name, ext_unit_dict): + if model.verbose: + print(f" loading {name} ...") + prop_array = [0] * nlay + for layer in range(nlay): + util2d_shape = get_util2d_shape_for_layer(model, layer=layer) + prop_array[layer] = Util2d.load( + f_obj, model, util2d_shape, dtype, name, ext_unit_dict + ) + return prop_array + + @staticmethod + def _ftype(): + return "DPF" + + @staticmethod + def _defaultunit(): + return 157 diff --git a/flopy/mfusg/mfusgdpt.py b/flopy/mfusg/mfusgdpt.py new file mode 100644 index 0000000000..c4185938eb --- /dev/null +++ b/flopy/mfusg/mfusgdpt.py @@ -0,0 +1,616 @@ +""" +Mfusgdpt module. + +Contains the MfUsgDpt class. Note that the user can +access the MfUsgDpt class as `flopy.mfusg.MfUsgDpt`. +""" + +import numpy as np + +from ..pakbase import Package +from ..utils import Util2d, Util3d +from ..utils.utils_def import ( + get_open_file_object, + get_unitnumber_from_ext_unit_dict, + get_util2d_shape_for_layer, +) +from .mfusg import MfUsg + + +class MfUsgDpt(Package): + """Dual Porosity Transport (dpt) Package Class for MODFLOW-USG Transport. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.modflow.Modflow`) to which + this package will be added. + ipakcb : int (0,1,-1), (default is 0) + a flag and a unit number >0 for cell-by-cell mass flux terms. + idptcon : int (0,1), (default is 0) + a flag and a unit number >0 for immobile domain concentrations + icbndimflg : int (0,1), (default is 1) + a flag active domain for the immobile (matrix) domain the same as + that for the mobile (fracture) domain + iadsorbim : int (0,1,2,3), (default is 0) + a flag for adsorption in the immobile domain(0: no adsorption, + 1: linear isotherm, 2: Freundlich isotherm, 3: Langmuir isotherm) + idispim : int (0,1), (default is 0) + a flag for dispersion in the immobile domain (0: no dispersion, + 1: dispersion) + izodim : int (0,1), (default is 0) + a flag for zero-order decay in the immobile domain (0: no zero-order decay, + 1: in water, 2: on soil, 3: on water and soil,4: on air-water interface) + ifodim : int (0,1), (default is 0) + a flag for first-order decay in the immobile domain (0: no first-order decay + 1: in water, 2: on soil, 3: on water and soil,4: on air-water interface) + frahk : bool, (default is False) + a flag for fractional hydraulic conductivity in the immobile domain + mobilesat : bool, (default is False) + immobile domain saturation equal to initial mobile domain saturation + inputsat : bool, (default is False) + a flag of immobile domain saturation input + icbundim : int or array of ints (nlay, nrow, ncol) + is cell-by-cell flag for transport simulation in immobile domain + phif : float or array of floats (nlay, nrow, ncol) + fraction of the total space that is occupied by the mobile domain + prsityim : float or array of floats (nlay, nrow, ncol) + effective transport porosity in the immobile domain + bulkdim : float or array of floats (nlay, nrow, ncol) + bulk density in the immobile domain + dlim : float or array of floats (nlay, nrow, ncol) + longitudinal dispersivity coefficient between mobile and immobile domains + ddtr : float or array of floats (nlay, nrow, ncol) + mass transfer coefficient between mobile and immobile domains + sim : float or array of floats (nlay, nrow, ncol) + saturation in the immobile domain + htcapsim : float or array of floats (nlay, nrow, ncol) + heat capacity in the immobile domain + htcondsim : float or array of floats (nlay, nrow, ncol) + heat conductivity in the immobile domain + adsorbim : float or array of floats (nlay, nrow, ncol) + adsorption coefficient in the immobile domain + flichim : float or array of floats (nlay, nrow, ncol) + Freundlich coefficient in the immobile domain + zodrwim : float or array of floats (nlay, nrow, ncol) + zero-order decay rate in water in the immobile domain + zodrsim : float or array of floats (nlay, nrow, ncol) + zero-order decay rate on soil in the immobile domain + fodrwim : float or array of floats (nlay, nrow, ncol) + first-order decay rate in water in the immobile domain + fodrsim : float or array of floats (nlay, nrow, ncol) + first-order decay rate on soil in the immobile domain + concim : float or array of floats (nlay, nrow, ncol) + initial concentration in the immobile domain + extension : string, (default is 'dpt'). + unitnumber : int, default is 58. + File unit number. + filenames : str or list of str + Filenames to use for the package and the output files. + add_package : bool, default is True + Flag to add the initialised package object to the parent model object. + + Methods + ------- + + See Also + -------- + + Notes + ----- + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> disu = flopy.mfusg.MfUsgDisU(model=ml, nlay=1, nodes=1, + iac=[1], njag=1,ja=np.array([0]), fahl=[1.0], cl12=[1.0]) + >>> dpt = flopy.mfusg.MfUsgdpt(ml)""" + + def __init__( + self, + model, + ipakcb=0, + idptcon=0, + icbndimflg=1, + iadsorbim=0, + idispim=0, + izodim=0, + ifodim=0, + frahk=False, + mobilesat=False, + inputsat=False, + icbundim=1, + phif=0.4, + prsityim=0.4, + bulkdim=1.6, + dlim=0.5, + ddtr=0.5, + sim=0.5, + htcapsim=0.0, + htcondsim=0.0, + adsorbim=0.0, + flichim=0.0, + zodrwim=0.0, + zodrsim=0.0, + fodrwim=0.0, + fodrsim=0.0, + concim=0.0, + extension="dpt", + unitnumber=None, + filenames=None, + add_package=True, + ): + """Constructs the MfUsgdpt object.""" + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + # set default unit number of one is not specified + if unitnumber is None: + self.unitnumber = self._defaultunit() + + super().__init__( + model, + extension, + self._ftype(), + unitnumber, + self._prepare_filenames(filenames), + ) + + self._generate_heading() + self.ipakcb = ipakcb + self.idptcon = idptcon + self.icbndimflg = icbndimflg + self.iadsorbim = iadsorbim + self.idispim = idispim + self.izodim = izodim + self.ifodim = ifodim + # options + self.frahk = frahk + self.mobilesat = mobilesat + self.inputsat = inputsat + + nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper + + if self.icbndimflg == 0: + self.icbundim = Util3d( + model, (nlay, nrow, ncol), np.int32, icbundim, name="icbundim" + ) + + self.phif = Util3d(model, (nlay, nrow, ncol), np.float32, phif, name="phif") + + self.prsityim = Util3d( + model, (nlay, nrow, ncol), np.float32, prsityim, name="prsityim" + ) + + if self.iadsorbim: + self.bulkdim = Util3d( + model, (nlay, nrow, ncol), np.float32, bulkdim, name="bulkdim" + ) + + self.dlim = Util3d(model, (nlay, nrow, ncol), np.float32, dlim, name="dlim") + + self.ddtr = Util3d(model, (nlay, nrow, ncol), np.float32, ddtr, name="ddtr") + + if self.inputsat: + self.sim = Util3d(model, (nlay, nrow, ncol), np.float32, sim, name="sim") + + if model.iheat: + self.htcapsim = Util3d( + model, (nlay, nrow, ncol), np.float32, htcapsim, name="htcapsim" + ) + + self.htcondsim = Util3d( + model, (nlay, nrow, ncol), np.float32, htcondsim, name="htcondsim" + ) + + mcomp = model.mcomp + + if isinstance(adsorbim, (int, float)): + adsorbim = [adsorbim] * mcomp + if isinstance(flichim, (int, float)): + flichim = [flichim] * mcomp + if isinstance(zodrwim, (int, float)): + zodrwim = [zodrwim] * mcomp + if isinstance(zodrsim, (int, float)): + zodrsim = [zodrsim] * mcomp + if isinstance(fodrwim, (int, float)): + fodrwim = [fodrwim] * mcomp + if isinstance(fodrsim, (int, float)): + fodrsim = [fodrsim] * mcomp + + if isinstance(concim, (int, float)): + concim = [concim] * mcomp + + self.adsorbim = [0] * mcomp + self.flichim = [0] * mcomp + self.zodrwim = [0] * mcomp + self.zodrsim = [0] * mcomp + self.fodrwim = [0] * mcomp + self.fodrsim = [0] * mcomp + self.concim = [0] * mcomp + + for icomp in range(mcomp): + if self.iadsorbim: + self.adsorbim[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + adsorbim[icomp], + name="adsorbim", + ) + + if self.iadsorbim == 2 or self.iadsorbim == 3: + self.flichim[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + flichim[icomp], + name="flichim", + ) + + if self.izodim == 1 or self.izodim == 3 or self.izodim == 4: + self.zodrwim[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + zodrwim[icomp], + name="zodrwim", + ) + + if self.iadsorbim and ( + self.izodim == 2 or self.izodim == 3 or self.izodim == 4 + ): + self.zodrsim[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + zodrsim[icomp], + name="zodrsim", + ) + + if self.ifodim == 1 or self.ifodim == 3 or self.ifodim == 4: + self.fodrwim[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + fodrwim[icomp], + name="fodrwim", + ) + + if self.iadsorbim and ( + self.ifodim == 2 or self.ifodim == 3 or self.ifodim == 4 + ): + self.fodrsim[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + fodrsim[icomp], + name="fodrsim", + ) + + self.concim[icomp] = Util3d( + model, (nlay, nrow, ncol), np.float32, concim[icomp], name="concim" + ) + + if add_package: + self.parent.add_package(self) + + def write_file(self, f=None): + """ + Write the dpt package file. + + Parameters + ---------- + f : open file object. + Default is None, which will result in MfUsg.fn_path being + opened for writing. + + Examples + -------- + """ + # Open file for writing + if f is None: + f_obj = open(self.fn_path, "w") + + # f_obj.write(f"{self.heading}\n") + + # Item 0: IPAKCB, IDPTCON + f_obj.write( + f" {self.ipakcb:9d} {self.idptcon:9d} {self.icbndimflg:9d}" + f" {self.iadsorbim:9d} {self.idispim:9d} {self.izodim:9d} {self.ifodim:9d}" + ) + + # Options + if self.frahk: + f_obj.write(" FRAHK") + + if self.mobilesat: + f_obj.write(" MOBILESAT") + + if self.inputsat: + f_obj.write(" INPUTSAT") + + f_obj.write("\n") + + # Item 1: ICBUNDIM + if self.icbndimflg == 0: + f_obj.write(self.icbundim.get_file_entry()) + + # Item 2: PHIF + f_obj.write(self.phif.get_file_entry()) + + # Item 3: PRSITYIM + f_obj.write(self.prsityim.get_file_entry()) + + # Item 4: BULKDIM + if self.iadsorbim: + f_obj.write(self.bulkdim.get_file_entry()) + + # Item 5: DLIM + if self.parent.idpf: + f_obj.write(self.dlim.get_file_entry()) + + # Item 6: DDTR + f_obj.write(self.ddtr.get_file_entry()) + + # Item 7: SIM + if self.inputsat: + f_obj.write(self.sim.get_file_entry()) + + # Item 8: HTCAPSIM, HTCONDSIM + if self.parent.iheat == 1: + f_obj.write(self.htcapsim.get_file_entry()) + f_obj.write(self.htcondsim.get_file_entry()) + + # Item 9: ADSORBIM, FLICHIM, ZODRWIM, ZODRSIM, FODRWIM, FODRSIM, CONCIM + mcomp = self.parent.mcomp + for icomp in range(mcomp): + if self.iadsorbim: + f_obj.write(self.adsorbim[icomp].get_file_entry()) + + if self.iadsorbim == 2 or self.iadsorbim == 3: + f_obj.write(self.flichim[icomp].get_file_entry()) + + if self.izodim == 1 or self.izodim == 3 or self.izodim == 4: + f_obj.write(self.zodrwim[icomp].get_file_entry()) + + if self.iadsorbim and ( + self.izodim == 2 or self.izodim == 3 or self.izodim == 4 + ): + f_obj.write(self.zodrsim[icomp].get_file_entry()) + + if self.ifodim == 1 or self.ifodim == 3 or self.ifodim == 4: + f_obj.write(self.fodrwim[icomp].get_file_entry()) + + if self.iadsorbim and ( + self.ifodim == 2 or self.ifodim == 3 or self.ifodim == 4 + ): + f_obj.write(self.fodrsim[icomp].get_file_entry()) + + f_obj.write(self.concim[icomp].get_file_entry()) + + # close the file + f_obj.close() + + @classmethod + def load(cls, f, model, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + dpt : MfUsgdpt object + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> dis = flopy.mfusg.MfUsgDisU.load('SeqDegEg.dis', ml) + >>> dpt = flopy.mfusg.MfUsgdpt.load('SeqDegEg.dpt', ml) + """ + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if model.verbose: + print("loading dpt package file...") + + f_obj = get_open_file_object(f, "r") + + # determine problem dimensions + nlay = model.nlay + + # item 0 + line = f_obj.readline().upper() + while line.startswith("#"): + line = f_obj.readline().upper() + + t = line.split() + kwargs = {} + + # item 1a + vars = { + "ipakcb": int, + "idptcon": int, + "icbndimflg": int, + "iadsorbim": int, + "idispim": int, + "izodim": int, + "ifodim": int, + } + + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = c(t[i].strip()) + # print(f"{v}={kwargs[v]}") + + # item 1a - options + if "frahk" in t: + kwargs["frahk"] = 1 + else: + kwargs["frahk"] = 0 + + if "mobilesat" in t: + kwargs["mobilesat"] = 1 + else: + kwargs["mobilesat"] = 0 + + if "inputsat" in t: + kwargs["inputsat"] = 1 + else: + kwargs["inputsat"] = 0 + + # item 1b + if kwargs["icbndimflg"] == 0: + kwargs["icbundim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.int32, "icbundim", ext_unit_dict + ) + + # item 2 + kwargs["phif"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "phif", ext_unit_dict + ) + + # item 3 + kwargs["prsityim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "prsityim", ext_unit_dict + ) + + # item 4 + if kwargs["iadsorbim"]: + kwargs["bulkdim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "bulkdim", ext_unit_dict + ) + + # item 5 + if model.idpf: + kwargs["dlim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "dlim", ext_unit_dict + ) + + # item 6 + kwargs["ddtr"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "ddtr", ext_unit_dict + ) + + # item 7 + if kwargs["inputsat"]: + kwargs["sim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "sim", ext_unit_dict + ) + + # item 8 + if model.iheat == 1: + kwargs["htcapsim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "htcapsim", ext_unit_dict + ) + + kwargs["htcondsim"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "htcondsim", ext_unit_dict + ) + + # item 9 + mcomp = model.mcomp + adsorbim = [0] * mcomp + flichim = [0] * mcomp + zodrwim = [0] * mcomp + zodrsim = [0] * mcomp + fodrwim = [0] * mcomp + fodrsim = [0] * mcomp + concim = [0] * mcomp + + for icomp in range(mcomp): + if kwargs["iadsorbim"]: + adsorbim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "adsorbim", ext_unit_dict + ) + + if kwargs["iadsorbim"] == 2 or kwargs["iadsorbim"] == 3: + flichim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "flichim", ext_unit_dict + ) + + if kwargs["izodim"] == 1 or kwargs["izodim"] == 3 or kwargs["izodim"] == 4: + zodrwim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "zodrwim", ext_unit_dict + ) + + if kwargs["iadsorbim"] and ( + kwargs["izodim"] == 2 or kwargs["izodim"] == 3 or kwargs["izodim"] == 4 + ): + zodrsim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "zodrsim", ext_unit_dict + ) + + if kwargs["ifodim"] == 1 or kwargs["ifodim"] == 3 or kwargs["ifodim"] == 4: + fodrwim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "fodrwim", ext_unit_dict + ) + + if kwargs["iadsorbim"] and ( + kwargs["ifodim"] == 2 or kwargs["ifodim"] == 3 or kwargs["ifodim"] == 4 + ): + fodrsim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "fodrsim", ext_unit_dict + ) + + concim[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "concim", ext_unit_dict + ) + + kwargs["adsorbim"] = adsorbim + kwargs["flichim"] = flichim + kwargs["zodrwim"] = zodrwim + kwargs["zodrsim"] = zodrsim + kwargs["fodrwim"] = fodrwim + kwargs["fodrsim"] = fodrsim + kwargs["concim"] = concim + + f_obj.close() + # set package unit number + unitnumber, filenames = get_unitnumber_from_ext_unit_dict( + model, cls, ext_unit_dict, kwargs["ipakcb"] + ) + + return cls(model, unitnumber=unitnumber, filenames=filenames, **kwargs) + + @staticmethod + def _load_prop_arrays(f_obj, model, nlay, dtype, name, ext_unit_dict): + if model.verbose: + print(f" loading {name} ...") + prop_array = [0] * nlay + for layer in range(nlay): + util2d_shape = get_util2d_shape_for_layer(model, layer=layer) + prop_array[layer] = Util2d.load( + f_obj, model, util2d_shape, dtype, name, ext_unit_dict + ) + return prop_array + + @staticmethod + def _ftype(): + return "DPT" + + @staticmethod + def _defaultunit(): + return 158 diff --git a/flopy/mfusg/mfusglak.py b/flopy/mfusg/mfusglak.py new file mode 100644 index 0000000000..532ddeccb0 --- /dev/null +++ b/flopy/mfusg/mfusglak.py @@ -0,0 +1,924 @@ +""" +mfusglak module. Contains the MfUsgLak class. Note that the user can access +the MfUsgLak class as `flopy.mfusg.MfUsgLak`. + +""" + +import numpy as np + +from ..pakbase import Package +from ..utils import Util3d, read_fixed_var, utils_def, write_fixed_var +from ..utils.util_array import Transient3d + + +class MfUsgLak(Package): + """ + MODFLOW USG Transport Lake Package Class. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.mfusg.MfUsg`) to which + this package will be added. + nlakes : int + NLAKES Number of separate lakes. + Sublakes of multiple-lake systems are considered separate lakes for + input purposes. The variable NLAKES is used, with certain internal + assumptions and approximations, to dimension arrays for the simulation. + ipakcb : int, optional + Toggles whether cell-by-cell budget data should be saved. If None or zero, + budget data will not be saved (default is None). + lwrt : int or list of ints (one per SP) + lwrt > 0, suppresses printout from the lake package. Default is 0 (to + print budget information) + theta : float + Explicit (THETA = 0.0), semi-implicit (0.0 < THETA < 1.0), or implicit + (THETA = 1.0) solution for lake stages. SURFDEPTH is read only if + THETA is assigned a negative value (the negative value of THETA is + then changed to a positive value internally by the code). + * A new method of solving for lake stage uses only the time-weighting + factor THETA (Merritt and Konikow, 2000, p. 52) for transient + simulations. THETA is automatically set to a value of 1.0 for all + steady-state stress periods. For transient stress periods, Explicit + (THETA = 0.0), semi-implicit (0.0 < THETA < 1.0), or implicit + (THETA = 1.0) solutions can be used to calculate lake stages. The + option to specify negative values for THETA is supported to allow + specification of additional variables (NSSITER, SSCNCR, SURFDEP) + for simulations that only include transient stress periods. If + THETA is specified as a negative value, then it is converted to a + positive value for calculations of lake stage. + * In MODFLOW-2000 and later, ISS is not part of the input. Instead + NSSITR or SSCNCR should be included if one or more stress periods + is a steady state stress period as defined in Ss/tr in the + Discretization file. + * SSCNCR and NSSITR can be read for a transient only simulation by + placing a negative sign immediately in front of THETA. A negative + THETA sets a flag which assumes input values for NSSITR and SSCNCR + will follow THETA in the format as described by Merritt and Konikow + (p. 52). A negative THETA is automatically reset to a positive + value after values of NSSITR and SSCNCR are read. + nssitr : int + Maximum number of iterations for Newton's method of solution for + equilibrium lake stages in each MODFLOW iteration for steady-state + aquifer head solution. Only read if ISS (option flag input to DIS + Package of MODFLOW indicating steady-state solution) is not zero or + if THETA is specified as a negative value. + * NSSITR and SSCNCR may be omitted for transient solutions (ISS = 0). + * In MODFLOW-2000 and later, ISS is not part of the input. + Instead NSSITR or SSCNCR should be included if one or more stress + periods is a steady state stress period as defined in Ss/tr in the + Discretization file. + * SSCNCR and NSSITR can be read for a transient only simulation by + placing a negative sign immediately in front of THETA. A negative + THETA sets a flag which assumes input values for NSSITR and SSCNCR + will follow THETA in the format as described by Merritt and Konikow + (p. 52). A negative THETA is automatically reset to a positive + value after values of NSSITR and SSCNCR are read. + * If NSSITR = 0, a value of 100 will be used instead. + sscncr : float + Convergence criterion for equilibrium lake stage solution by Newton's + method. Only read if ISS is not zero or if THETA is specified as a + negative value. See notes above for nssitr. + surfdepth : float + The height of small topological variations (undulations) in lake-bottom + elevations that can affect groundwater discharge to lakes. SURFDEPTH + decreases the lakebed conductance for vertical flow across a horizontal + lakebed caused both by a groundwater head that is between the lakebed + and the lakebed plus SURFDEPTH and a lake stage that is also between + the lakebed and the lakebed plus SURFDEPTH. This method provides a + smooth transition from a condition of no groundwater discharge to a + lake, when groundwater head is below the lakebed, to a condition of + increasing groundwater discharge to a lake as groundwater head becomes + greater than the elevation of the dry lakebed. The method also allows + for the transition of seepage from a lake to groundwater when the lake + stage decreases to the lakebed elevation. Values of SURFDEPTH ranging + from 0.01 to 0.5 have been used successfully in test simulations. + SURFDEP is read only if THETA is specified as a negative value. + stages : float or list of floats + The initial stage of each lake at the beginning of the run. + stage_range : list of tuples (ssmn, ssmx) of length nlakes + Where ssmn and ssmx are the minimum and maximum stages allowed for each + lake in steady-state solution. + * SSMN and SSMX are not needed for a transient run and must be + omitted when the solution is transient. + * When the first stress period is a steady-state stress period, + SSMN is defined in record 3. + + For subsequent steady-state stress periods, SSMN is defined in + record 9a. + lakarr : array of integers (nlay, nrow, ncol) + LKARR A value is read in for every grid cell. + If LKARR(I,J,K) = 0, the grid cell is not a lake volume cell. + If LKARR(I,J,K) > 0, its value is the identification number of the lake + occupying the grid cell. LKARR(I,J,K) must not exceed the value NLAKES. + If it does, or if LKARR(I,J,K) < 0, LKARR(I,J,K) is set to zero. + Lake cells cannot be overlain by non-lake cells in a higher layer. + Lake cells must be inactive cells (IBOUND = 0) and should not be + convertible to active cells (WETDRY = 0). + + The Lake package can be used when all or some of the model layers + containing the lake are confined. The authors recommend using the + Layer-Property Flow Package (LPF) for this case, although the + BCF and HUF Packages will work too. However, when using the BCF6 + package to define aquifer properties, lake/aquifer conductances in the + lateral direction are based solely on the lakebed leakance (and not on + the lateral transmissivity of the aquifer layer). As before, when the + BCF6 package is used, vertical lake/aquifer conductances are based on + lakebed conductance and on the vertical hydraulic conductivity of the + aquifer layer underlying the lake when the wet/dry option is + implemented, and only on the lakebed leakance when the wet/dry option + is not implemented. + bdlknc : array of floats (nlay, nrow, ncol) + BDLKNC A value is read in for every grid cell. The value is the lakebed + leakance that will be assigned to lake/aquifer interfaces that occur + in the corresponding grid cell. If the wet-dry option flag (IWDFLG) is + not active (cells cannot rewet if they become dry), then the BDLKNC + values are assumed to represent the combined leakances of the lakebed + material and the aquifer material between the lake and the centers of + the underlying grid cells, i. e., the vertical conductance values (CV) + will not be used in the computation of conductances across lake/aquifer + boundary faces in the vertical direction. + + IBOUND and WETDRY should be set to zero for every cell for which LKARR + is not equal to zero. IBOUND is defined in the input to the Basic + Package of MODFLOW. WETDRY is defined in the input to the BCF or other + flow package of MODFLOW if the IWDFLG option is active. When used with + the HUF package, the Lake Package has been modified to compute + effective lake-aquifer conductance solely on the basis of the + user-specified value of lakebed leakance; aquifer hydraulic + conductivities are not used in this calculation. An appropriate + informational message is now printed after the lakebed conductances + are written to the main output file. + sill_data : dict + (dataset 8 in documentation) + Dict of lists keyed by stress period. Each list has a tuple of dataset + 8a, 8b for every multi-lake system, where dataset 8a is another tuple of + IC : int + The number of sublakes + ISUB : list of ints + The identification numbers of the sublakes in the sublake + system being described in this record. The center lake number + is listed first. + And dataset 8b contains + SILLVT : sequence of floats + A sequence of sill elevations for each sublakes that determines + whether the center lake is connected with a given sublake. + Values are entered for each sublake in the order the sublakes + are listed in the previous record. + flux_data : dict + (dataset 9 in documentation) + Dict of lists keyed by stress period. The list for each stress period + is a list of lists, with each list containing the variables + PRCPLK EVAPLK RNF WTHDRW [SSMN] [SSMX] from the documentation. + PRCPLK : float + The rate of precipitation per unit area at the surface of a + lake (L/T). + EVAPLK : float + The rate of evaporation per unit area from the surface of a + lake (L/T). + RNF : float + Overland runoff from an adjacent watershed entering the lake. + If RNF > 0, it is specified directly as a volumetric rate, or + flux (L3 /T). If RNF < 0, its absolute value is used as a + dimensionless multiplier applied to the product of the lake + precipitation rate per unit area (PRCPLK) and the surface area + of the lake at its full stage (occupying all layer 1 lake + cells). When RNF is entered as a dimensionless multiplier + (RNF < 0), it is considered to be the product of two + proportionality factors. The first is the ratio of the area of + the basin contributing runoff to the surface area of the lake + when it is at full stage. The second is the fraction of the + current rainfall rate that becomes runoff to the lake. This + procedure provides a means for the automated computation of + runoff rate from a watershed to a lake as a function of + varying rainfall rate. For example, if the basin area is 10 + times greater than the surface area of the lake, and 20 percent + of the precipitation on the basin becomes overland runoff + directly into the lake, then set RNF = -2.0. + WTHDRW : float + The volumetric rate, or flux (L3 /T), of water removal from a + lake by means other than rainfall, evaporation, surface + outflow, or groundwater seepage. A negative value indicates + augmentation. Normally, this would be used to specify the + rate of artificial withdrawal from a lake for human water use, + or if negative, artificial augmentation of a lake volume for + aesthetic or recreational purposes. + SSMN : float + Minimum stage allowed for each lake in steady-state solution. + See notes on ssmn and ssmx above. + SSMX : float + SSMX Maximum stage allowed for each lake in steady-state + solution. + + options : list of strings + Package options. (default is None). + extension : string + Filename extension (default is 'lak') + unitnumber : int + File unit number (default is None). + filenames : str or list of str + Filenames to use for the package and the output files. If + filenames=None the package name will be created using the model name + and package extension and the cbc output name will be created using + the model name and .cbc extension (for example, modflowtest.cbc), + if ipakcb is a number greater than zero. If a single string is passed + the package will be set to the string and cbc output names will be + created using the model name and .cbc extension, if ipakcb is a + number greater than zero. To define the names for all package files + (input and output) the length of the list of strings should be 2. + Default is None. + + Methods + ------- + + See Also + -------- + + Notes + ----- + Parameters are not supported in FloPy. + + Examples + -------- + + >>> import flopy + >>> m = flopy.mfusg.MfUsg() + >>> lak = {} + >>> lak[0] = [[2, 3, 4, 15.6, 1050., -4]] #this lake boundary will be + >>> #applied to all stress periods + >>> lak = flopy.mfusg.MfUsgLak(m, nstress_period_data=strd) + + """ + + def __init__( + self, + model, + nlakes=1, + ipakcb=None, + theta=-1.0, + nssitr=0, + sscncr=0.0, + surfdep=0.0, + stages=1.0, + stage_range=None, + clake=None, + tab_files=None, + tab_units=None, + lakarr=None, + bdlknc=None, + sill_data=None, + flux_data=None, + conc_data=None, + extension="lak", + unitnumber=None, + filenames=None, + options=None, + lwrt=0, + **kwargs, + ): + # set default unit number of one is not specified + if unitnumber is None: + unitnumber = MfUsgLak._defaultunit() + + # set filenames + tabdata = False + nlen = 2 + if options is not None: + for option in options: + if "TABLEINPUT" in option.upper(): + tabdata = True + nlen += nlakes + break + filenames = self._prepare_filenames(filenames, nlen) + + # cbc output file + self.set_cbc_output_file(ipakcb, model, filenames[1]) + + # table input files + if tabdata: + if tab_files is None: + tab_files = filenames[2:] + + # add tab_files as external files + if tabdata: + # make sure the number of tabfiles is equal to the number of lakes + if len(tab_files) < nlakes: + msg = ( + "a tabfile must be specified for each lake " + "{} tabfiles specified " + "instead of {} tabfiles".format(len(tab_files), nlakes) + ) + # TODO: what should happen with msg? + # make sure tab_files are not None + for idx, fname in enumerate(tab_files, 1): + if fname is None: + raise ValueError( + "a filename must be specified for the " + f"tabfile for lake {idx}" + ) + # set unit for tab files if not passed to __init__ + if tab_units is None: + tab_units = [] + for idx in range(len(tab_files)): + tab_units.append(model.next_ext_unit()) + # add tabfiles as external files + for iu, fname in zip(tab_units, tab_files): + model.add_external(fname, iu) + + # call base package constructor + super().__init__( + model, + extension=extension, + name=self._ftype(), + unit_number=unitnumber, + filenames=filenames[0], + ) + + self._generate_heading() + self.url = "lak.html" + + if options is None: + options = [] + self.options = options + self.nlakes = nlakes + self.theta = theta + self.nssitr = nssitr + self.sscncr = sscncr + self.surfdep = surfdep + self.lwrt = lwrt + + if isinstance(stages, float): + if self.nlakes == 1: + stages = np.array([self.nlakes], dtype=float) * stages + else: + stages = np.ones(self.nlakes, dtype=float) * stages + elif isinstance(stages, list): + stages = np.array(stages) + if stages.shape[0] != nlakes: + err = f"stages shape should be ({nlakes}) but is only ({stages.shape[0]})." + raise Exception(err) + self.stages = stages + if stage_range is None: + stage_range = np.ones((nlakes, 2), dtype=float) + stage_range[:, 0] = -10000.0 + stage_range[:, 1] = 10000.0 + else: + if isinstance(stage_range, list): + stage_range = np.array(stage_range) + elif isinstance(stage_range, float): + raise Exception( + f"stage_range should be a list or array of size ({nlakes}, 2)" + ) + + self.dis = utils_def.get_dis(model) + if self.dis.steady[0]: + if stage_range.shape != (nlakes, 2): + raise Exception( + "stages shape should be ({},2) but is only {}.".format( + nlakes, stage_range.shape + ) + ) + self.stage_range = stage_range + + # tabfile data + self.tabdata = tabdata + self.iunit_tab = tab_units + + if lakarr is None and bdlknc is None: + err = "lakarr and bdlknc must be specified" + raise Exception(err) + nrow, ncol, nlay, nper = self.parent.get_nrow_ncol_nlay_nper() + self.lakarr = Transient3d( + model, (nlay, nrow, ncol), np.int32, lakarr, name="lakarr_" + ) + self.bdlknc = Transient3d( + model, (nlay, nrow, ncol), np.float32, bdlknc, name="bdlknc_" + ) + + if sill_data is not None: + if not isinstance(sill_data, dict): + try: + sill_data = {0: sill_data} + except: + err = "sill_data must be a dictionary" + raise Exception(err) + + if flux_data is not None: + if not isinstance(flux_data, dict): + # convert array to a dictionary + try: + flux_data = {0: flux_data} + except: + err = "flux_data must be a dictionary" + raise Exception(err) + for key, value in flux_data.items(): + if isinstance(value, np.ndarray): + td = {} + for k in range(value.shape[0]): + td[k] = value[k, :].tolist() + flux_data[key] = td + if len(list(flux_data.keys())) != nlakes: + raise Exception( + f"flux_data dictionary must have {nlakes} entries" + ) + elif isinstance(value, float) or isinstance(value, int): + td = {} + for k in range(self.nlakes): + td[k] = (np.ones(6, dtype=float) * value).tolist() + flux_data[key] = td + elif isinstance(value, dict): + try: + steady = self.dis.steady[key] + except: + steady = True + nlen = 4 + if steady and key > 0: + nlen = 6 + for k in range(self.nlakes): + td = value[k] + if len(td) < nlen: + raise Exception( + "flux_data entry for stress period {} " + "has {} entries but should have " + "{} entries".format(key + 1, nlen, len(td)) + ) + + self.flux_data = flux_data + self.sill_data = sill_data + + self.clake = clake + self.conc_data = conc_data + + self.parent.add_package(self) + + return + + def _ncells(self): + """Maximum number of cells that can have lakes (developed for + MT3DMS SSM package). + + Returns + ------- + ncells: int + maximum number of lak cells + + """ + nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper + return nlay * nrow * ncol + + def write_file(self): + """ + Write the package file. + + Returns + ------- + None + + """ + f = open(self.fn_path, "w") + # dataset 0 + if self.parent.version != "mf2k": + f.write(f"{self.heading}\n") + + # dataset 1a + if len(self.options) > 0: + for option in self.options: + f.write(f"{option} ") + f.write("\n") + + # dataset 1b + f.write( + write_fixed_var( + [self.nlakes, self.ipakcb], free=self.parent.free_format_input + ) + ) + # dataset 2 + steady = np.any(self.dis.steady.array) + t = [self.theta] + if self.theta < 0.0 or steady: + t.append(self.nssitr) + t.append(self.sscncr) + if self.theta < 0.0: + t.append(self.surfdep) + f.write(write_fixed_var(t, free=self.parent.free_format_input)) + + # dataset 3 + mcomp = self.parent.mcomp + + steady = self.dis.steady[0] + for n in range(self.nlakes): + ipos = [10] + t = [self.stages[n]] + if steady: + ipos.append(10) + t.append(self.stage_range[n, 0]) + ipos.append(10) + t.append(self.stage_range[n, 1]) + if mcomp > 0: + for icomp in range(mcomp): + ipos.append(10) + t.append(self.clake[n][icomp]) + if self.tabdata: + ipos.append(5) + t.append(self.iunit_tab[n]) + f.write(write_fixed_var(t, ipos=ipos, free=self.parent.free_format_input)) + + ds8_keys = list(self.sill_data.keys()) if self.sill_data is not None else [] + ds9_keys = list(self.flux_data.keys()) + nper = self.dis.steady.shape[0] + for kper in range(nper): + itmp, file_entry_lakarr = self.lakarr.get_kper_entry(kper) + ibd, file_entry_bdlknc = self.bdlknc.get_kper_entry(kper) + + itmp2 = 0 + if kper in ds9_keys: + itmp2 = 1 + elif len(ds9_keys) > 0: + itmp2 = -1 + if isinstance(self.lwrt, list): + tmplwrt = self.lwrt[kper] + else: + tmplwrt = self.lwrt + t = [itmp, itmp2, tmplwrt] + comment = f"Stress period {kper + 1}" + f.write( + write_fixed_var(t, free=self.parent.free_format_input, comment=comment) + ) + + if itmp > 0: + f.write(file_entry_lakarr) + f.write(file_entry_bdlknc) + + nslms = 0 + if kper in ds8_keys: + ds8 = self.sill_data[kper] + nslms = len(ds8) + + f.write( + write_fixed_var( + [nslms], + length=5, + free=self.parent.free_format_input, + comment="Data set 7", + ) + ) + if nslms > 0: + for n in range(nslms): + d1, d2 = ds8[n] + s = write_fixed_var( + d1, + length=5, + free=self.parent.free_format_input, + comment="Data set 8a", + ) + f.write(s) + s = write_fixed_var( + d2, + free=self.parent.free_format_input, + comment="Data set 8b", + ) + f.write(s) + + if itmp2 > 0: + ds9 = self.flux_data[kper] + ds9b = self.conc_data[kper] + for n in range(self.nlakes): + try: + steady = self.dis.steady[kper] + except: + steady = True + if kper > 0 and steady: + t = ds9[n] + else: + t = ds9[n][0:4] + s = write_fixed_var( + t, + free=self.parent.free_format_input, + comment="Data set 9a", + ) + f.write(s) + if mcomp > 0: + for icomp in range(mcomp): + t = ds9b[n, icomp] + s = write_fixed_var( + t, + free=self.parent.free_format_input, + comment="Data set 9b", + ) + f.write(s) + + # close the lak file + f.close() + + @classmethod + def load(cls, f, model, nper=None, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + nper : int + The number of stress periods. If nper is None, then nper will be + obtained from the model object. (default is None). + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + lak : ModflowLak object + ModflowLak object. + + Examples + -------- + + >>> import flopy + >>> m = flopy.mfusg.MfUsg() + >>> lak = flopy.mfusg.MfUsgLak.load('test.lak', m) + + """ + + cls.dis = utils_def.get_dis(model) + + if model.verbose: + print("loading lak package file...") + + openfile = not hasattr(f, "read") + if openfile: + filename = f + f = open(filename, "r", errors="replace") + + # dataset 0 -- header + while True: + line = f.readline() + if line[0] != "#": + break + + options = [] + tabdata = False + if "TABLEINPUT" in line.upper(): + if model.verbose: + print(" reading lak dataset 1a") + options.append("TABLEINPUT") + tabdata = True + line = f.readline() + + # read dataset 1b + if model.verbose: + print(" reading lak dataset 1b") + t = line.strip().split() + nlakes = int(t[0]) + ipakcb = 0 + try: + ipakcb = int(t[1]) + except: + pass + + # read dataset 2 + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + t = read_fixed_var(line, ncol=4) + theta = float(t[0]) + nssitr, sscncr = 0, 0.0 + if theta < 0: + try: + nssitr = int(t[1]) + except: + if model.verbose: + print(" implicit nssitr defined in file") + try: + sscncr = float(t[2]) + except: + if model.verbose: + print(" implicit sscncr defined in file") + + surfdep = 0.0 + if theta < 0.0: + surfdep = float(t[3]) + + if nper is None: + nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper() + + mcomp = model.mcomp + + if model.verbose: + print(" reading lak dataset 3") + stages = [] + stage_range = [] + clake = [] + if tabdata: + tab_units = [] + else: + tab_units = None + for lake in range(nlakes): + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + t = read_fixed_var(line, ipos=[10, 10, 10, 5]) + stages.append(t[0]) + ipos = 1 + if cls.dis.steady[0]: + stage_range.append((float(t[ipos]), float(t[ipos + 1]))) + ipos += 2 + if mcomp > 0: + conc = [] + for icomp in range(mcomp): + conc.append(float(t[ipos])) + ipos += 1 + clake.append(conc) + if tabdata: + iu = int(t[ipos]) + tab_units.append(iu) + + lake_loc = {} + lake_lknc = {} + sill_data = {} + flux_data = {} + conc_data = {} + lwrt = [] + for iper in range(nper): + if model.verbose: + print(f" reading lak dataset 4 - for stress period {iper + 1}") + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + t = read_fixed_var(line, ncol=3) + itmp, itmp1, tmplwrt = int(t[0]), int(t[1]), int(t[2]) + lwrt.append(tmplwrt) + + if itmp > 0: + if model.verbose: + print(f" reading lak dataset 5 - for stress period {iper + 1}") + name = f"LKARR_StressPeriod_{iper}" + lakarr = Util3d.load( + f, model, (nlay, nrow, ncol), np.int32, name, ext_unit_dict + ) + if model.verbose: + print(f" reading lak dataset 6 - for stress period {iper + 1}") + name = f"BDLKNC_StressPeriod_{iper}" + bdlknc = Util3d.load( + f, + model, + (nlay, nrow, ncol), + np.float32, + name, + ext_unit_dict, + ) + + lake_loc[iper] = lakarr + lake_lknc[iper] = bdlknc + + if model.verbose: + print(f" reading lak dataset 7 - for stress period {iper + 1}") + line = f.readline().rstrip() + t = line.split() + nslms = int(t[0]) + ds8 = [] + if nslms > 0: + if model.verbose: + print( + f" reading lak dataset 8 - for stress period {iper + 1}" + ) + for i in range(nslms): + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + ic = int(line[0:5]) + t = read_fixed_var(line, ncol=ic + 1, length=5) + ic = int(t[0]) + ds8a = [ic] + for j in range(1, ic + 1): + ds8a.append(int(t[j])) + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + t = read_fixed_var(line, ncol=ic - 1) + silvt = [] + for j in range(ic - 1): + silvt.append(float(t[j])) + ds8.append((ds8a, silvt)) + sill_data[iper] = ds8 + if itmp1 >= 0: + if model.verbose: + print(f" reading lak dataset 9 - for stress period {iper + 1}") + ds9 = {} + ds9b = {} + for n in range(nlakes): + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + t = read_fixed_var(line, ncol=6) + tds = [] + tds.append(float(t[0])) # PRCPLK + tds.append(float(t[1])) # EVAPLK + tds.append(float(t[2])) # RNF + tds.append(float(t[3])) # WTHDRW + + wthdrw = float(t[3]) + + if cls.dis.steady[iper]: + if iper == 0: + tds.append(stage_range[n][0]) + tds.append(stage_range[n][1]) + else: + tds.append(float(t[4])) + tds.append(float(t[5])) + else: + tds.append(0.0) + tds.append(0.0) + ds9[n] = tds + + if mcomp > 0: + for icomp in range(mcomp): + line = f.readline().rstrip() + if model.array_free_format: + t = line.split() + else: + t = read_fixed_var(line, ncol=3) + tds = [] + tds.append(float(t[0])) # CPPT + tds.append(float(t[1])) # CRNF + if wthdrw < 0: + tds.append(float(t[2])) # CAUG + ds9b[(n, icomp)] = tds + flux_data[iper] = ds9 + conc_data[iper] = ds9b + + if openfile: + f.close() + + # convert lake data to Transient3d objects + lake_loc = Transient3d( + model, (nlay, nrow, ncol), np.int32, lake_loc, name="lakarr_" + ) + lake_lknc = Transient3d( + model, (nlay, nrow, ncol), np.float32, lake_lknc, name="bdlknc_" + ) + + # determine specified unit number + n = 2 + if tab_units is not None: + n += nlakes + unitnumber = None + filenames = [None for x in range(n)] + if ext_unit_dict is not None: + unitnumber, filenames[0] = model.get_ext_dict_attr( + ext_unit_dict, filetype=cls._ftype() + ) + if ipakcb > 0: + iu, filenames[1] = model.get_ext_dict_attr(ext_unit_dict, unit=ipakcb) + model.add_pop_key_list(ipakcb) + + ipos = 2 + if tab_units is not None: + for i in range(len(tab_units)): + iu, filenames[ipos] = model.get_ext_dict_attr( + ext_unit_dict, unit=tab_units[i] + ) + ipos += 1 + + return cls( + model, + options=options, + nlakes=nlakes, + ipakcb=ipakcb, + theta=theta, + nssitr=nssitr, + surfdep=surfdep, + sscncr=sscncr, + lwrt=lwrt, + stages=stages, + stage_range=stage_range, + clake=clake, + tab_units=tab_units, + lakarr=lake_loc, + bdlknc=lake_lknc, + sill_data=sill_data, + flux_data=flux_data, + conc_data=conc_data, + unitnumber=unitnumber, + filenames=filenames, + ) + + @staticmethod + def _ftype(): + return "LAK" + + @staticmethod + def _defaultunit(): + return 119 diff --git a/flopy/mfusg/mfusgmdt.py b/flopy/mfusg/mfusgmdt.py new file mode 100644 index 0000000000..5b22f5a6e3 --- /dev/null +++ b/flopy/mfusg/mfusgmdt.py @@ -0,0 +1,454 @@ +""" +Mfusgmdt module. + +Contains the MfUsgMdt class. Note that the user can +access the MfUsgMdt class as `flopy.mfusg.MfUsgMdt`. +""" + +import numpy as np + +from ..pakbase import Package +from ..utils import Util2d, Util3d +from ..utils.utils_def import ( + get_open_file_object, + get_unitnumber_from_ext_unit_dict, + get_util2d_shape_for_layer, +) +from .mfusg import MfUsg + + +class MfUsgMdt(Package): + """Matrix Diffusion Transport (mdt) Package Class for MODFLOW-USG Transport. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.modflow.Modflow`) to which + this package will be added. + ipakcb : int (0,1,-1), (default is 0) + a flag and a unit number >0 for cell-by-cell mass flux terms. + imdtcf : int (0,1,2,3,4,5,6,7), (default is 0) + a flag and a unit number >0 for immobile domain concentrations + mdflag : int, (default is 0) + a flag for the matrix diffusion type + volfracmd : float, (default is 0.0) + volume fraction of fracture/high permeability domain. + SAME AS PHIF IN DPT PACKAGE + pormd : float, (default is 0.0) + effective transport porosity of the mobile domain + rhobmd : float, (default is 0.0) + bulk density of the mobile domain + difflenmd : float, (default is 0.0) + diffusion length for diffusive transport within the matrix domain + tortmd : float, (default is 0.0) + tortuosity factor for the matrix domain + kdmd : float, (default is 0.0) + adsorption coefficient (Kd) of a contaminant species in immobile domain + decaymd : float, (default is 0.0) + first-order decay rate of a contaminant species in the immobile domain + yieldmd : float, (default is 0.0) + yield coefficient for chain decay + diffmd : float, (default is 0.0) + diffusion coefficient for the immobile domain + aiold1md : float, (default is 0.0) (MDFLAG = 1, 3, 4, 5, 6, or 7) + concentration integral (equation 8) of each species + aiold2md : float, (default is 0.0) (MDFLAG = 2, 5, 6, or 7) + concentration integral (equation 8) of each species + frahk : bool, (default is False) + True = hydraulic conductivity and storage terms only for mobile domain + False = hydraulic conductivity and storage terms for matrix domains + fradarcy : bool, (default is False) + True = Darcy flux computed only for fracture domain + tshiftmd : float, (default is 0.0) + time shift for the immobile domain + iunitAI2 :float, (default is 0.0) + separate binary file of ai2 for the matrix domain + crootname_md : str, (default is None) + rootname for the concentration output binary file + extension : str, (default is 'mdt'). + unitnumber : int, (default is 57). + File unit number and the output files. + filenames : str or list of str + Filenames to use for the package and the output files. + add_package : bool, default is True + Flag to add the initialised package object to the parent model object. + + Methods + ------- + + See Also + -------- + + Notes + ----- + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg(exe_name='USGs_1.exe') + >>> disu = flopy.mfusg.MfUsgDisU(model=ml, nlay=1, nodes=1, + iac=[1], njag=1,ja=np.array([0]), fahl=[1.0], cl12=[1.0]) + >>> mdt = flopy.mfusg.MfUsgmdt(ml)""" + + def __init__( + self, + model, + ipakcb=0, + imdtcf=0, + mdflag=0, + volfracmd=0.0, + pormd=0.0, + rhobmd=0.0, + difflenmd=0.0, + tortmd=0.0, + kdmd=0.0, + decaymd=0.0, + yieldmd=0.0, + diffmd=0.0, + aiold1md=0.0, + aiold2md=0.0, + frahk=False, + fradarcy=False, + tshiftmd=0.0, + iunitAI2=0, + crootname=None, + extension="mdt", + unitnumber=None, + filenames=None, + add_package=True, + ): + """Constructs the MfUsgmdt object.""" + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + # set default unit number of one is not specified + if unitnumber is None: + self.unitnumber = self._defaultunit() + + super().__init__( + model, + extension, + self._ftype(), + unitnumber, + self._prepare_filenames(filenames), + ) + + self._generate_heading() + self.ipakcb = ipakcb + self.imdtcf = imdtcf + + # options + self.frahk = frahk + self.fradarcy = fradarcy + self.tshiftmd = tshiftmd + self.iunitAI2 = iunitAI2 + self.crootname = crootname + + nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper + + self.mdflag = Util3d(model, (nlay, nrow, ncol), np.int32, mdflag, name="mdflag") + self.volfracmd = Util3d( + model, (nlay, nrow, ncol), np.float32, volfracmd, name="volfracmd" + ) + self.pormd = Util3d(model, (nlay, nrow, ncol), np.float32, pormd, name="pormd") + self.rhobmd = Util3d( + model, (nlay, nrow, ncol), np.float32, rhobmd, name="rhobmd" + ) + self.difflenmd = Util3d( + model, (nlay, nrow, ncol), np.float32, difflenmd, name="difflenmd" + ) + self.tortmd = Util3d( + model, (nlay, nrow, ncol), np.float32, tortmd, name="tortmd" + ) + + mcomp = model.mcomp + + self.kdmd = [0] * mcomp + self.decaymd = [0] * mcomp + self.yieldmd = [0] * mcomp + self.diffmd = [0] * mcomp + self.aiold1md = [0] * mcomp + self.aiold2md = [0] * mcomp + + for icomp in range(mcomp): + self.kdmd[icomp] = Util3d( + model, (nlay, nrow, ncol), np.float32, kdmd[icomp], name="kdmd" + ) + self.decaymd[icomp] = Util3d( + model, (nlay, nrow, ncol), np.float32, decaymd[icomp], name="decaymd" + ) + self.yieldmd[icomp] = Util3d( + model, (nlay, nrow, ncol), np.float32, yieldmd[icomp], name="yieldmd" + ) + self.diffmd[icomp] = Util3d( + model, (nlay, nrow, ncol), np.float32, diffmd[icomp], name="diffmd" + ) + if self.tshiftmd > 0: + self.aiold1md[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + aiold1md[icomp], + name="aiold1md", + ) + self.aiold2md[icomp] = Util3d( + model, + (nlay, nrow, ncol), + np.float32, + aiold2md[icomp], + name="aiold2md", + ) + + if add_package: + self.parent.add_package(self) + + def write_file(self, f=None): + """ + Write the mdt package file. + + Parameters + ---------- + f : open file object. + Default is None, which will result in MfUsg.fn_path being + opened for writing. + + Examples + -------- + """ + # Open file for writing + if f is None: + f_obj = open(self.fn_path, "w") + + f_obj.write(f"{self.heading}\n") + f_obj.write(f"{ self.ipakcb:9d} {self.imdtcf:9d}") + + # Options + if self.frahk: + f_obj.write(" FRAHK") + + if self.fradarcy: + f_obj.write(" FRADARCY") + + if self.tshiftmd > 0: + f_obj.write(f" TSHIFTMD {self.tshiftmd:9.2f}") + + if self.iunitAI2 > 0: + f_obj.write(f" SEPARATE_AI2 {self.iunitAI2:9d}") + + if self.crootname is not None: + f_obj.write(f" MULTIFILE_MD {self.crootname}") + + f_obj.write("\n") + + # Item 1: MDFLAG, VOLFRACMD, PORMD, RHOBMD, DIFFLENMD, TORTMD + f_obj.write(self.mdflag.get_file_entry()) + if not self.parent.idpf: + f_obj.write(self.volfracmd.get_file_entry()) + f_obj.write(self.pormd.get_file_entry()) + f_obj.write(self.rhobmd.get_file_entry()) + f_obj.write(self.difflenmd.get_file_entry()) + f_obj.write(self.tortmd.get_file_entry()) + + # Item 2 - 7: KDM, DECAYMD, YIELDMD, DIFFMD, AIOLD1MD, AIOLD2MD + mcomp = self.parent.mcomp + + for icomp in range(mcomp): + f_obj.write(self.kdmd[icomp].get_file_entry()) + f_obj.write(self.decaymd[icomp].get_file_entry()) + f_obj.write(self.yieldmd[icomp].get_file_entry()) + f_obj.write(self.diffmd[icomp].get_file_entry()) + if self.tshiftmd > 0: + f_obj.write(self.aiold1md[icomp].get_file_entry()) + f_obj.write(self.aiold2md[icomp].get_file_entry()) + + # close the file + f_obj.close() + + @classmethod + def load(cls, f, model, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + mdt : MfUsgmdt object + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> dis = flopy.mfusg.MfUsgDisU.load('SeqDegEg.dis', ml) + >>> mdt = flopy.mfusg.MfUsgmdt.load('SeqDegEg.mdt', ml) + """ + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if model.verbose: + print("loading mdt package file...") + + f_obj = get_open_file_object(f, "r") + + nlay = model.nlay + + # item 0 + line = f_obj.readline().upper() + while line[0] == "#": + line = f_obj.readline().upper() + + print(f"line={line}") + + t = line.split() + kwargs = {} + + # item 1a + vars = {"ipakcb": int, "imdtcf": int} + + for i, (v, c) in enumerate(vars.items()): + kwargs[v] = c(t[i].strip()) + # print(f"{v}={kwargs[v]}\n") + + # item 1a - options + if "frahk" in t: + kwargs["frahk"] = 1 + else: + kwargs["frahk"] = 0 + + if "fradarcy" in t: + kwargs["fradarcy"] = 1 + else: + kwargs["fradarcy"] = 0 + + if "TSHIFTMD" in t: + idx = t.index("TSHIFTMD") + kwargs["tshiftmd"] = float(t[idx + 1]) + else: + kwargs["tshiftmd"] = 0.0 + + if "SEPARATE_AI2" in t: + idx = t.index("SEPARATE_AI2") + kwargs["iunitAI2"] = int(t[idx + 1]) + else: + kwargs["iunitAI2"] = 0 + + if "MULTIFILE_MD" in t: + idx = t.index("MULTIFILE_MD") + kwargs["crootname"] = str(t[idx + 1]) + else: + kwargs["crootname"] = None + + kwargs["mdflag"] = cls._load_prop_arrays( + f_obj, model, nlay, np.int32, "mdflag", ext_unit_dict + ) + + if not model.idpf: + kwargs["volfracmd"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "volfracmd", ext_unit_dict + ) + + kwargs["pormd"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "pormd", ext_unit_dict + ) + + kwargs["rhobmd"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "rhobmd", ext_unit_dict + ) + + kwargs["difflenmd"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "difflenmd", ext_unit_dict + ) + + kwargs["tortmd"] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "tortmd", ext_unit_dict + ) + + mcomp = model.mcomp + + kdmd = [0] * mcomp + decaymd = [0] * mcomp + yieldmd = [0] * mcomp + diffmd = [0] * mcomp + aiold1md = [0] * mcomp + aiold2md = [0] * mcomp + + for icomp in range(mcomp): + kdmd[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "kdmd", ext_unit_dict + ) + + decaymd[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "decaymd", ext_unit_dict + ) + + yieldmd[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "yieldmd", ext_unit_dict + ) + + diffmd[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "diffmd", ext_unit_dict + ) + + if kwargs["tshiftmd"] > 0: + aiold1md[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "aiold1md", ext_unit_dict + ) + + aiold2md[icomp] = cls._load_prop_arrays( + f_obj, model, nlay, np.float32, "aiold2md", ext_unit_dict + ) + + kwargs["kdmd"] = kdmd + kwargs["decaymd"] = decaymd + kwargs["yieldmd"] = yieldmd + kwargs["diffmd"] = diffmd + kwargs["aiold1md"] = aiold1md + kwargs["aiold2md"] = aiold2md + + f_obj.close() + # set package unit number + unitnumber, filenames = get_unitnumber_from_ext_unit_dict( + model, cls, ext_unit_dict, kwargs["ipakcb"] + ) + + return cls(model, unitnumber=unitnumber, filenames=filenames, **kwargs) + + @staticmethod + def _load_prop_arrays(f_obj, model, nlay, dtype, name, ext_unit_dict): + if model.verbose: + print(f" loading {name} ...") + prop_array = [0] * nlay + for layer in range(nlay): + util2d_shape = get_util2d_shape_for_layer(model, layer=layer) + prop_array[layer] = Util2d.load( + f_obj, model, util2d_shape, dtype, name, ext_unit_dict + ) + return prop_array + + @staticmethod + def _ftype(): + return "MDT" + + @staticmethod + def _defaultunit(): + return 156 diff --git a/flopy/mfusg/mfusgoc.py b/flopy/mfusg/mfusgoc.py new file mode 100644 index 0000000000..865e27a842 --- /dev/null +++ b/flopy/mfusg/mfusgoc.py @@ -0,0 +1,1240 @@ +""" +mfusgoc module. Contains the MfusgOc class. Note that the user can access +the MfusgOc class as `flopy.mfusg.MfusgOc`. + +""" + +import os + +from ..pakbase import Package +from ..utils.utils_def import type_from_iterable +from .mfusg import MfUsg + + +class MfUsgOc(Package): + """ + MODFLOW USG Transport Output Control (OC) Package Class. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.mfusg.MfUsg`) to which + this package will be added. + ihedfm : int (default is 0) + is a code for the format in which heads will be printed. + iddnfm : int (default is 0) + is a code for the format in which drawdown will be printed. + chedfm : string (default is None) + is a character value that specifies the format for saving heads. + cddnfm : string (default is None) + is a character value that specifies the format for saving drawdown. + cboufm : string (default is None) + is a character value that specifies the format for saving ibound. + stress_period_data : dictionary of lists + Dictionary key is a tuple with the zero-based period and step + (IPEROC, ITSOC) for each print/save option list. If stress_period_data + is None, then heads are saved for the last time step of each stress + period. (default is None) + + The list can have any valid MFUSG OC print/save option: + PRINT HEAD + PRINT DRAWDOWN + PRINT BUDGET + SAVE HEAD + SAVE DRAWDOWN + SAVE BUDGET + SAVE IBOUND + + The lists can also include (1) DDREFERENCE in the list to reset + drawdown reference to the period and step and (2) a list of layers + for PRINT HEAD, SAVE HEAD, PRINT DRAWDOWN, SAVE DRAWDOWN, and + SAVE IBOUND. + + stress_period_data = {(0,1):['save head']}) would save the head for + the second timestep in the first stress period. + + compact : boolean + Save results in compact budget form. (default is True). + extension : list of strings + (default is ['oc', 'hds', 'ddn', 'cbc', 'ibo']). + unitnumber : list of ints + (default is [14, 51, 52, 53, 0]). + filenames : str or list of str + + Attributes + ---------- + + Methods + ------- + + See Also + -------- + + Notes + ----- + The "words" method for specifying output control is the only option + available. Also, the "compact" budget should normally be used as it + produces files that are typically much smaller. The compact budget form is + also a requirement for using the MODPATH particle tracking program. + + Examples + -------- + + >>> import flopy + >>> m = flopy.mfusg.MfUsg() + >>> spd = {(0, 0): ['print head'], + ... (0, 1): [], + ... (0, 249): ['print head'], + ... (0, 250): [], + ... (0, 499): ['print head', 'save ibound'], + ... (0, 500): [], + ... (0, 749): ['print head', 'ddreference'], + ... (0, 750): [], + ... (0, 999): ['print head']} + >>> oc = flopy.mfusg.MfusgOc(m, stress_period_data=spd, cboufm='(20i5)') + + """ + + def __init__( + self, + model, + ihedfm=0, + iddnfm=0, + ispcfm=0, + chedfm=None, + cddnfm=None, + cspcfm=None, + cboufm=None, + compact=True, + stress_period_data={(0, 0): ["save head"]}, + extension=["oc", "hds", "ddn", "cbc", "ibo", "con"], + unitnumber=None, + filenames=None, + label="LABEL", + **kwargs, + ): + """Constructs the MfusgOc object.""" + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if unitnumber is None: + unitnumber = MfUsgOc._defaultunit() + elif isinstance(unitnumber, list): + if len(unitnumber) < 6: + for idx in range(len(unitnumber), 6): + unitnumber.append(0) + self.label = label + + # set filenames + filenames = self._prepare_filenames(filenames, 6) + + # support structured and unstructured dis + dis = model.get_package("DIS") + if dis is None: + dis = model.get_package("DISU") + + if stress_period_data is None: + stress_period_data = { + (kper, dis.nstp.array[kper] - 1): ["save head"] + for kper in range(dis.nper) + } + + # process kwargs + if "save_every" in kwargs: + save_every = int(kwargs.pop("save_every")) + else: + save_every = None + if save_every is not None: + if "save_types" in kwargs: + save_types = kwargs.pop("save_types") + if isinstance(save_types, str): + save_types = [save_types] + else: + save_types = ["save head", "print budget"] + if "save_conc" in kwargs: + save_types.append("save conc") + if "save_start" in kwargs: + save_start = int(kwargs.pop("save_start")) + else: + save_start = 1 + stress_period_data = {} + for kper in range(dis.nper): + icnt = save_start + for kstp in range(dis.nstp[kper]): + if icnt == save_every: + stress_period_data[(kper, kstp)] = save_types + icnt = 0 + else: + stress_period_data[(kper, kstp)] = [] + icnt += 1 + + # adaptive time stepping + self.atsa = 0 + self.nptimes = 0 + self.npsteps = 0 + self.timot = None + + # FASTFORWARD initial conditions of heads + self.fastforward = False + self.ispfast = 0 + self.itsfast = 0 + self.iugfast = 0 + self.iucfast = 0 + self.iudfast = 0 + + # FASTFORWARDC initial conditions of concentration + self.fastforwardc = False + self.ispfastc = 0 + self.itsfastc = 0 + self.iugfastc = 0 + self.iucfastc = 0 + self.iudfastc = 0 + self.iumfastc = 0 + + if "atsa" in kwargs: + self.atsa = int(kwargs.pop("atsa")) + if "nptimes" in kwargs: + self.nptimes = int(kwargs.pop("nptimes")) + if "npsteps" in kwargs: + self.npsteps = int(kwargs.pop("npsteps")) + if "timot" in kwargs: + self.timot = kwargs.pop("timot") + + if "fastforward" in kwargs: + self.fastforward = int(kwargs.pop("fastforward")) + + if self.fastforward: + if "ispfast" in kwargs: + self.ispfast = int(kwargs.pop("ispfast")) + if "itsfast" in kwargs: + self.itsfast = int(kwargs.pop("itsfast")) + if "iugfast" in kwargs: + self.iugfast = int(kwargs.pop("iugfast")) + if "iucfast" in kwargs: + self.iucfast = int(kwargs.pop("iucfast")) + if "iudfast" in kwargs: + self.iudfast = int(kwargs.pop("iudfast")) + + if "fastforwardc" in kwargs: + self.fastforwardc = int(kwargs.pop("fastforwardc")) + + if self.fastforwardc: + if "ispfastc" in kwargs: + self.ispfastc = int(kwargs.pop("ispfastc")) + if "itsfastc" in kwargs: + self.itsfastc = int(kwargs.pop("itsfastc")) + if "iugfastc" in kwargs: + self.iugfastc = int(kwargs.pop("iugfastc")) + if "iucfastc" in kwargs: + self.iucfastc = int(kwargs.pop("iucfastc")) + if "iudfastc" in kwargs: + self.iudfastc = int(kwargs.pop("iudfastc")) + if "iumfastc" in kwargs: + self.iumfastc = int(kwargs.pop("iumfastc")) + + # set output unit numbers based on oc settings + self.savehead = False + self.saveddn = False + self.savebud = False + self.saveibnd = False + self.savespc = False + for key, value in stress_period_data.items(): + for t in list(value): + tlwr = t.lower() + if "save head" in tlwr: + self.savehead = True + if unitnumber[1] == 0: + unitnumber[1] = 51 + if "save drawdown" in tlwr: + self.saveddn = True + if unitnumber[2] == 0: + unitnumber[2] = 52 + if "save budget" in tlwr: + self.savebud = True + if unitnumber[3] == 0 and filenames is None: + unitnumber[3] = 53 + if "save ibound" in tlwr: + self.saveibnd = True + if unitnumber[4] == 0: + unitnumber[4] = 54 + if "save conc" in tlwr: + self.savespc = True + if unitnumber[5] == 0: + unitnumber[5] = 55 + + # do not create head, ddn, or cbc output files if output is not + # specified in the oc stress_period_data + if not self.savehead: + unitnumber[1] = 0 + if not self.saveddn: + unitnumber[2] = 0 + if not self.savebud: + unitnumber[3] = 0 + if not self.saveibnd: + unitnumber[4] = 0 + if not self.savespc: + unitnumber[5] = 0 + + self.iuhead = unitnumber[1] + self.iuddn = unitnumber[2] + self.iubud = unitnumber[3] + self.iuibnd = unitnumber[4] + self.iuspc = unitnumber[5] + + # add output files + # head file + if self.savehead: + model.add_output_file( + unitnumber[1], + fname=filenames[1], + extension=extension[1], + binflag=chedfm is None, + ) + # drawdown file + if self.saveddn: + model.add_output_file( + unitnumber[2], + fname=filenames[2], + extension=extension[2], + binflag=cddnfm is None, + ) + # budget file + # Nothing is needed for the budget file + + # ibound file + if self.saveibnd: + model.add_output_file( + unitnumber[4], + fname=filenames[4], + extension=extension[4], + binflag=cboufm is None, + ) + + # concentration file + if self.savespc: + model.add_output_file( + unitnumber[5], + fname=filenames[5], + extension=extension[5], + binflag=cspcfm is None, + ) + + # call base package constructor + super().__init__( + model, + extension=extension[0], + name=self._ftype(), + unit_number=unitnumber[0], + filenames=filenames[0], + ) + + self._generate_heading() + + self.url = "oc.html" + self.ihedfm = ihedfm + self.iddnfm = iddnfm + self.chedfm = chedfm + self.cddnfm = cddnfm + + self.cboufm = cboufm + + self.ispcfm = ispcfm + self.cspcfm = cspcfm + + self.compact = compact + + self.stress_period_data = stress_period_data + + self.parent.add_package(self) + + def check(self, f=None, verbose=True, level=1, checktype=None): + """ + Check package data for common errors. + + Parameters + ---------- + f : str or file handle + String defining file name or file handle for summary file + of check method output. If a string is passed a file handle + is created. If f is None, check method does not write + results to a summary file. (default is None) + verbose : bool + Boolean flag used to determine if check method results are + written to the screen. + level : int + Check method analysis level. If level=0, summary checks are + performed. If level=1, full checks are performed. + + Returns + ------- + None + + Examples + -------- + + >>> import flopy + >>> m = flopy.mfusg.MfUsg.load('model.nam') + >>> m.oc.check() + + """ + chk = self._get_check(f, verbose, level, checktype) + dis = self.parent.get_package("DIS") + if dis is None: + dis = self.parent.get_package("DISU") + if dis is None: + chk._add_to_summary("Error", package="OC", desc="DIS package not available") + else: + # generate possible actions expected + expected_actions = [] + for first in ["PRINT", "SAVE"]: + for second in ["HEAD", "DRAWDOWN", "BUDGET", "IBOUND", "CONC"]: + expected_actions.append([first, second]) + # remove exception + del expected_actions[expected_actions.index(["PRINT", "IBOUND"])] + keys = list(self.stress_period_data.keys()) + for kper in range(dis.nper): + for kstp in range(dis.nstp[kper]): + kperkstp = (kper, kstp) + if kperkstp in keys: + del keys[keys.index(kperkstp)] + data = self.stress_period_data[kperkstp] + if not isinstance(data, list): + data = [data] + for action in data: + words = action.upper().split() + if len(words) < 2: + chk._add_to_summary( + "Warning", + package="OC", + desc=f"action {action!r} ignored; too few words", + ) + elif words[0:2] not in expected_actions: + chk._add_to_summary( + "Warning", + package="OC", + desc=f"action {action!r} ignored", + ) + # TODO: check data list of layers for some actions + for kperkstp in keys: + # repeat as many times as remaining keys not used + chk._add_to_summary( + "Warning", + package="OC", + desc="action(s) defined in OC stress_period_data ignored " + "as they are not part the stress periods defined by DIS", + ) + chk.summarize() + return chk + + def write_file(self): + """ + Write the package file. + + Returns + ------- + None + + """ + f_oc = open(self.fn_path, "w") + f_oc.write(f"{self.heading}\n") + + # write options + if self.atsa: + f_oc.write("ATSA ") + if self.nptimes > 0: + f_oc.write(f"NPTIMES {self.nptimes} ") + for i in range(self.nptimes): + f_oc.write(f"{self.timot[i]} ") + if self.npsteps > 0: + f_oc.write(f"NPSTPS {self.npsteps} ") + f_oc.write("\n") + if self.fastforward: + f_oc.write( + f"FASTFORWARD {self.ispfast:3.0f} {self.itsfast:3.0f}" + f" {self.iugfast:3.0f} {self.iucfast:3.0f} {self.iudfast:3.0f}\n" + ) + if self.fastforwardc: + f_oc.write( + f"FASTFORWARDC {self.ispfastc:3.0f} {self.itsfastc:3.0f}" + f" {self.iugfastc:3.0f} {self.iucfastc:3.0f} {self.iudfastc:3.0f}" + f" {self.iumfastc:3.0f}\n" + ) + + if self.atsa and self.nptimes > 0: + for i in range(self.nptimes): + f_oc.write(f"{self.timot[i]} ") + f_oc.write("\n") + + line = f"HEAD PRINT FORMAT {self.ihedfm:3.0f}\n" + f_oc.write(line) + if self.chedfm is not None: + line = f"HEAD SAVE FORMAT {self.chedfm:20s} {self.label}\n" + f_oc.write(line) + if self.savehead: + line = f"HEAD SAVE UNIT {self.iuhead:5.0f}\n" + f_oc.write(line) + + f_oc.write(f"DRAWDOWN PRINT FORMAT {self.iddnfm:3.0f}\n") + if self.cddnfm is not None: + line = f"DRAWDOWN SAVE FORMAT {self.cddnfm:20s} {self.label}\n" + f_oc.write(line) + if self.saveddn: + line = f"DRAWDOWN SAVE UNIT {self.iuddn:5.0f}\n" + f_oc.write(line) + + if self.saveibnd: + if self.cboufm is not None: + line = f"IBOUND SAVE FORMAT {self.cboufm:20s} {self.label}\n" + f_oc.write(line) + line = f"IBOUND SAVE UNIT {self.iuibnd:5.0f}\n" + f_oc.write(line) + + if self.savespc: + f_oc.write(f"CONC PRINT FORMAT {self.ispcfm:3.0f}\n") + if self.cspcfm is not None: + line = f"CONC SAVE FORMAT {self.cspcfm:20s} {self.label}\n" + f_oc.write(line) + line = f"CONC SAVE UNIT {self.iuspc:5.0f}\n" + f_oc.write(line) + + if self.compact: + f_oc.write("COMPACT BUDGET AUX\n") + + # add a line separator between header and stress + # period data + # f_oc.write("\n") + + # write the transient sequence described by the data dict + nr, nc, nl, nper = self.parent.get_nrow_ncol_nlay_nper() + dis = self.parent.get_package("DIS") + if dis is None: + dis = self.parent.get_package("DISU") + nstp = dis.nstp + + keys = list(self.stress_period_data.keys()) + keys.sort() + + data = [] + ddnref = "" + lines = "" + for kper in range(nper): + for kstp in range(nstp[kper]): + kperkstp = (kper, kstp) + if kperkstp in keys: + data = self.stress_period_data[kperkstp] + if not isinstance(data, list): + data = [data] + lines = "" + if len(data) > 0: + for item in data: + if "DDREFERENCE" in item.upper(): + ddnref = item.lower() + else: + lines += f" {item}\n" + if len(lines) > 0: + f_oc.write(f"period {kper + 1} step {kstp + 1} {ddnref}\n") + f_oc.write(lines) + # f_oc.write("\n") + ddnref = "" + lines = "" + + # close oc file + f_oc.close() + + def _set_singlebudgetunit(self, budgetunit): + if budgetunit is None: + budgetunit = self.parent.next_ext_unit() + self.iubud = budgetunit + + def _set_budgetunit(self): + iubud = [] + for i, pp in enumerate(self.parent.packagelist): + if hasattr(pp, "ipakcb"): + if pp.ipakcb > 0: + iubud.append(pp.ipakcb) + if len(iubud) < 1: + iubud = None + elif len(iubud) == 1: + iubud = iubud[0] + self.iubud = iubud + + def get_budgetunit(self): + """ + Get the budget file unit number(s). + + Parameters + ---------- + None + + Returns + ------- + iubud : integer or list of integers + Unit number or list of cell-by-cell budget output unit numbers. + None is returned if ipakcb is less than one for all packages. + """ + + # set iubud by iterating through the packages + self._set_budgetunit() + return self.iubud + + def reset_budgetunit(self, budgetunit=None, fname=None): + """ + Reset the cell-by-cell budget unit (ipakcb) for every package that + can write cell-by-cell data when SAVE BUDGET is specified in the + OC file to the specified budgetunit. + + Parameters + ---------- + budgetunit : int, optional + Unit number for cell-by-cell output data. If budgetunit is None + then the next available external unit number is assigned. Default + is None + fname : string, optional + Filename to use for cell-by-cell output file. If fname=None the + cell-by-cell output file will be created using the model name and + a '.cbc' file extension. Default is None. + + Returns + ------- + None + + """ + + # remove existing output file + for pp in self.parent.packagelist: + if hasattr(pp, "ipakcb"): + if pp.ipakcb > 0: + self.parent.remove_output(unit=pp.ipakcb) + pp.ipakcb = 0 + + # set the unit number used for all cell-by-cell output + self._set_singlebudgetunit(budgetunit) + + # add output file + for pp in self.parent.packagelist: + if hasattr(pp, "ipakcb"): + pp.ipakcb = self.iubud + self.parent.add_output_file(pp.ipakcb, fname=fname, package=pp.name) + + return + + @staticmethod + def get_ocoutput_units(f, ext_unit_dict=None): + """ + Get head and drawdown units from a OC file. + + Parameters + ---------- + f : filename or file handle + File to load. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + ihedun : integer + Unit number of the head file. + fhead : str + File name of the head file. Is only defined if ext_unit_dict is + passed and the unit number is a valid key. + , headfilename, oc : MfUsgOc object + MfUsgOc object. + iddnun : integer + Unit number of the drawdown file. + fddn : str + File name of the drawdown file. Is only defined if ext_unit_dict is + passed and the unit number is a valid key. + + Examples + -------- + + >>> import flopy + >>> ihds, hf, iddn, df = flopy.mfusg.MfUsgOc.get_ocoutput_units('test.oc') + + """ + + # initialize + ihedun = 0 + iddnun = 0 + fhead = None + fddn = None + + numericformat = False + + openfile = not hasattr(f, "read") + if openfile: + filename = f + f = open(filename, "r") + + # read header + ipos = f.tell() + while True: + line = f.readline() + if line[0] == "#": + continue + elif line[0] == []: + continue + else: + lnlst = line.strip().split() + try: + ihedun, iddnun = int(lnlst[2]), int(lnlst[3]) + numericformat = True + except: + f.seek(ipos) + # exit so the remaining data can be read + # from the file based on numericformat + break + # read word formats + if not numericformat: + while True: + line = f.readline() + if len(line) < 1: + break + lnlst = line.strip().split() + if line[0] == "#": + continue + + # skip blank line in the OC file + if len(lnlst) < 1: + continue + + # dataset 1 values + elif ( + "HEAD" in lnlst[0].upper() + and "SAVE" in lnlst[1].upper() + and "UNIT" in lnlst[2].upper() + ): + ihedun = int(lnlst[3]) + elif ( + "DRAWDOWN" in lnlst[0].upper() + and "SAVE" in lnlst[1].upper() + and "UNIT" in lnlst[2].upper() + ): + iddnun = int(lnlst[3]) + # dataset 2 + elif "PERIOD" in lnlst[0].upper(): + break + # + if ext_unit_dict is not None: + if ihedun in ext_unit_dict: + fhead = ext_unit_dict[ihedun] + if iddnun in ext_unit_dict: + fddn = ext_unit_dict[iddnun] + + if openfile: + f.close() + + # return + return ihedun, fhead, iddnun, fddn + + @classmethod + def load(cls, f, model, nper=None, nstp=None, nlay=None, ext_unit_dict=None): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.mfusg.MfUsg`) to + which this package will be added. + nper : int + The number of stress periods. If nper is None, then nper will be + obtained from the model object. (default is None). + nstp : int or list of ints + Integer of list of integers containing the number of time steps + in each stress period. If nstp is None, then nstp will be obtained + from the DIS or DISU packages attached to the model object. The + length of nstp must be equal to nper. (default is None). + nlay : int + The number of model layers. If nlay is None, then nnlay will be + obtained from the model object. nlay only needs to be specified + if an empty model object is passed in and the oc file being loaded + is defined using numeric codes. (default is None). + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + oc : MfUsgOc object + MfUsgOc object. + + Examples + -------- + + >>> import flopy + >>> m = flopy.mfusg.Mfusg() + >>> oc = flopy.mfusg.MfUsgOc.load('test.oc', m) + + """ + + if model.verbose: + print("loading oc package file...") + + # set nper + if nper is None or nlay is None: + nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper() + + if nper == 0 or nlay == 0: + raise ValueError( + "discretization package not defined for the model, " + "nper and nlay must be provided to the .load() method" + ) + + # set nstp + if nstp is None: + dis = model.get_package("DIS") + if dis is None: + dis = model.get_package("DISU") + if dis is None: + raise ValueError( + "discretization package not defined for the model, " + "a nstp list must be provided to the .load() method" + ) + nstp = list(dis.nstp.array) + else: + if isinstance(nstp, (int, float)): + nstp = [int(nstp)] + + # validate the size of nstp + if len(nstp) != nper: + raise OSError( + f"nstp must be a list with {nper} entries, " + f"provided nstp list has {len(nstp)} entries." + ) + + # initialize + ihedfm = 0 + iddnfm = 0 + + ihedun = 0 + iddnun = 0 + ibouun = 0 + + compact = False + chedfm = None + cddnfm = None + cboufm = None + + # concentration of species + ispcfm = 0 + ispcun = 0 + cspcfm = None + + stress_period_data = {} + + # read header + kwargs = {} + kwargs["atsa"] = 0 + kwargs["nptimes"] = 0 + kwargs["npsteps"] = 0 + kwargs["timot"] = None + + # FASTFORWARD initial conditions of heads + # stress period number to fast-forward + kwargs["ispfast"] = 0 + # time step number to fast-forward + kwargs["itsfast"] = 0 + # unit number for groundwater head file + kwargs["iugfast"] = 0 + # unit number for CLN domain head file + kwargs["iucfast"] = 0 + # unit number for dual domain head file + kwargs["iudfast"] = 0 + # FASTFORWARD initial conditions of concentration + # stress period number to fast-forward concentration + kwargs["ispfastc"] = 0 + # time step number to fast-forward concentration + kwargs["itsfastc"] = 0 + # unit number for groundwater concentration file + kwargs["iugfastc"] = 0 + # unit number for CLN domain concentration file + kwargs["iucfastc"] = 0 + # unit number for dual domain concentration file + kwargs["iudfastc"] = 0 + # unit number for matrix diffusion concentration file + kwargs["iumfastc"] = 0 + + # BOOTSTRAPPING initial conditions of heads + # unit number for groundwater head file + kwargs["iugboot"] = 0 + # unit number for CLN domain head file + kwargs["iucboot"] = 0 + # unit number for dual domain head file + kwargs["iudboot"] = 0 + + numericformat = False + + openfile = not hasattr(f, "read") + if openfile: + filename = f + f = open(filename, "r") + else: + filename = os.path.basename(f.name) + + line = f.readline().upper() + while line[0] == "#": + line = f.readline().upper() + + lnlst = line.strip().split() + try: + ihedfm, iddnfm = int(lnlst[0]), int(lnlst[1]) + ihedun, iddnun = int(lnlst[2]), int(lnlst[3]) + numericformat = True + except: + f.seek(0) + + # read options + ispcfm = type_from_iterable(lnlst, index=4) + ispcun = type_from_iterable(lnlst, index=5) + if "ATS" in lnlst or "ATSA" in lnlst: + kwargs["atsa"] = 1 + if "NPTIMES" in lnlst: + idx = lnlst.index("NPTIMES") + kwargs["nptimes"] = int(lnlst[idx + 1]) + if "NPSTPS" in lnlst: + idx = lnlst.index("NPSTPS") + kwargs["npsteps"] = int(lnlst[idx + 1]) + if "FASTFORWARD" in lnlst: + idx = lnlst.index("FASTFORWARD") + kwargs["ispfast"] = int(lnlst[idx + 1]) + kwargs["itsfast"] = int(lnlst[idx + 2]) + kwargs["iugfast"] = int(lnlst[idx + 3]) + if model.icln: + kwargs["iucfast"] = int(lnlst[idx + 4]) + if model.idpf: + kwargs["iudfast"] = int(lnlst[idx + 5]) + if "FASTFORWARDC" in lnlst: + idx = lnlst.index("FASTFORWARDC") + kwargs["ispfastc"] = int(lnlst[idx + 1]) + kwargs["itsfastc"] = int(lnlst[idx + 2]) + kwargs["iugfastc"] = int(lnlst[idx + 3]) + if model.icln: + kwargs["iucfastc"] = int(lnlst[idx + 4]) + if model.idpt: + kwargs["iudfastc"] = int(lnlst[idx + 5]) + if model.imdt: + kwargs["iumfastc"] = int(lnlst[idx + 5]) + if "BOOTSTRAPPING" in lnlst: + idx = lnlst.index("BOOTSTRAPPING") + kwargs["iugboot"] = int(lnlst[idx + 1]) + if model.icln: + kwargs["iucboot"] = int(lnlst[idx + 2]) + if model.idpf: + kwargs["iudboot"] = int(lnlst[idx + 3]) + + # read numeric formats + if numericformat: + # read TIMOT values + if kwargs["atsa"] and kwargs["nptimes"] > 0: + line = f.readline() + lnlst = line.strip().split() + kwargs["timot"] = [0] * kwargs["nptimes"] + for i in range(kwargs["nptimes"]): + kwargs["timot"][i] = int(lnlst[i]) + # process each line + lines = [] + # Item 3 is read for each time step if adaptive time stepping is not used. + for iperoc in range(nper): + # adaptive time stepping is used, read item 3 for each stress period + if kwargs["atsa"]: + nstp[iperoc] = 1 + for itsoc in range(nstp[iperoc]): + lines = [] + if kwargs["atsa"]: + line = f.readline() + lnlst = line.strip().split() + lines.append("DELTAT {float(lnlst[0]):11.4e}") + lines.append("TMINAT {float(lnlst[1]):11.4e}") + lines.append("TMAXAT {float(lnlst[2]):11.4e}") + lines.append("TADJAT {float(lnlst[3]):11.4e}") + lines.append("TCUTAT {float(lnlst[4]):11.4e}") + + line = f.readline() + lnlst = line.strip().split() + incode, ihddfl = int(lnlst[0]), int(lnlst[1]) + ibudfl, icbcfl = int(lnlst[2]), int(lnlst[3]) + ispcfl = type_from_iterable(lnlst, index=4) + + # new print and save flags are needed if incode is not + # less than 0. + if incode < 0: + if len(lines) > 0: + stress_period_data[(iperoc, itsoc)] = list(lines) + continue + # set print and save budget flags + if ibudfl != 0: + lines.append("PRINT BUDGET") + if icbcfl != 0: + lines.append("SAVE BUDGET") + if incode == 0: + line = f.readline() + lnlst = line.strip().split() + hdpr, ddpr = int(lnlst[0]), int(lnlst[1]) + hdsv, ddsv = int(lnlst[2]), int(lnlst[3]) + cnpr = type_from_iterable(lnlst, index=4) + cnsv = type_from_iterable(lnlst, index=5) + if hdpr != 0: + lines.append("PRINT HEAD") + if ddpr != 0: + lines.append("PRINT DRAWDOWN") + if cnpr != 0: + lines.append("PRINT CONC") + if hdsv != 0: + lines.append("SAVE HEAD") + if ddsv != 0: + lines.append("SAVE DRAWDOWN") + if cnsv != 0: + lines.append("SAVE CONC") + elif incode > 0: + headprint = "" + headsave = "" + ddnprint = "" + ddnsave = "" + spcprint = "" + spcsave = "" + for k in range(nlay): + line = f.readline() + lnlst = line.strip().split() + hdpr, ddpr = int(lnlst[0]), int(lnlst[1]) + hdsv, ddsv = int(lnlst[2]), int(lnlst[3]) + cnpr = type_from_iterable(lnlst, index=4) + cnsv = type_from_iterable(lnlst, index=5) + if hdpr != 0: + headprint += f" {k + 1}" + if ddpr != 0: + ddnprint += f" {k + 1}" + if hdsv != 0: + headsave += f" {k + 1}" + if ddsv != 0: + ddnsave += f" {k + 1}" + if cnpr != 0: + spcprint += f" {k + 1}" + if cnsv != 0: + spcsave += f" {k + 1}" + if len(headprint) > 0: + lines.append(f"PRINT HEAD{headprint}") + if len(ddnprint) > 0: + lines.append(f"PRINT DRAWDOWN{ddnprint}") + if len(headsave) > 0: + lines.append(f"SAVE HEAD{headsave}") + if len(ddnsave) > 0: + lines.append(f"SAVE DRAWDOWN{ddnsave}") + if len(spcprint) > 0: + lines.append(f"PRINT CONC{spcprint}") + if len(spcsave) > 0: + lines.append(f"SAVE CONC{spcsave}") + stress_period_data[(iperoc, itsoc)] = list(lines) + + # Output Control Using Words + else: + iperoc, itsoc = 0, 0 + lines = [] + while True: + line = f.readline().upper() + if not line: + if len(lines) > 0: + stress_period_data[(iperoc, itsoc)] = list(lines) + break + lnlst = line.strip().split() + if line[0] == "#": + continue + + # added by JJS 12/12/14 to avoid error when there is a blank line + if lnlst == []: + continue + # end add + + # dataset 1 values + + # TIMOT : output time values + elif lnlst[0].isdigit(): + if kwargs["atsa"] and kwargs["nptimes"] > 0: + kwargs["timot"] = [0] * kwargs["nptimes"] + for i in range(kwargs["nptimes"]): + kwargs["timot"][i] = int(lnlst[i]) + elif ( + "HEAD" in lnlst[0] and "PRINT" in lnlst[1] and "FORMAT" in lnlst[2] + ): + ihedfm = int(lnlst[3]) + + elif "HEAD" in lnlst[0] and "SAVE" in lnlst[1] and "FORMAT" in lnlst[2]: + chedfm = lnlst[3] + elif "HEAD" in lnlst[0] and "SAVE" in lnlst[1] and "UNIT" in lnlst[2]: + ihedun = int(lnlst[3]) + elif ( + "DRAWDOWN" in lnlst[0] + and "PRINT" in lnlst[1] + and "FORMAT" in lnlst[2] + ): + iddnfm = int(lnlst[3]) + elif ( + "DRAWDOWN" in lnlst[0] + and "SAVE" in lnlst[1] + and "FORMAT" in lnlst[2] + ): + cddnfm = lnlst[3] + elif ( + "DRAWDOWN" in lnlst[0] and "SAVE" in lnlst[1] and "UNIT" in lnlst[2] + ): + iddnun = int(lnlst[3]) + elif ( + "IBOUND" in lnlst[0] and "SAVE" in lnlst[1] and "FORMAT" in lnlst[2] + ): + cboufm = lnlst[3] + elif "IBOUND" in lnlst[0] and "SAVE" in lnlst[1] and "UNIT" in lnlst[2]: + ibouun = int(lnlst[3]) + elif ( + "CONC" in lnlst[0] and "PRINT" in lnlst[1] and "FORMAT" in lnlst[2] + ): + ispcfm = int(lnlst[3]) + elif "CONC" in lnlst[0] and "SAVE" in lnlst[1] and "FORMAT" in lnlst[2]: + cspcfm = lnlst[3] + elif "CONC" in lnlst[0] and "SAVE" in lnlst[1] and "UNIT" in lnlst[2]: + ispcun = int(lnlst[3]) + elif "COMPACT" in lnlst[0]: + compact = True + + # dataset 2 + elif "PERIOD" in lnlst[0]: + iperoc1 = int(lnlst[1]) - 1 + if "STEP" in lnlst: + itsoc1 = int(lnlst[3]) - 1 + else: + itsoc1 = 0 + + if len(lines) > 0: + if iperoc1 > 0 or itsoc1 > 0: + # create period step tuple + kperkstp = (iperoc, itsoc) + # save data + stress_period_data[kperkstp] = lines + # reset lines + lines = [] + + # update iperoc and itsoc + iperoc, itsoc = iperoc1, itsoc1 + + # do not used data that exceeds nper + if iperoc > nper: + break + + # dataset 3 + elif "PRINT" in lnlst[0]: + lines.append(f"{lnlst[0]} {lnlst[1]}") + elif "SAVE" in lnlst[0]: + lines.append(f"{lnlst[0]} {lnlst[1]}") + elif "DELTAT" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "TMINAT" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "TMAXAT" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "TADJAT" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "TCUTAT" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "HCLOSE" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "BTOL" in lnlst[0]: + lines.append(f"{lnlst[0]} {float(lnlst[1]):11.4e}") + elif "MXITER" in lnlst[0]: + lines.append(f"{lnlst[0]} {int(lnlst[1]):11d}") + elif "BOOTSTRAP" in lnlst[0]: + lines.append(f"{lnlst[0]}") + elif "NOBOOTSTRAP" in lnlst[0]: + lines.append(f"{lnlst[0]}") + elif "BOOTSTRAPSCALE" in lnlst[0]: + lines.append(f"{lnlst[0]}") + elif "NOBOOTSTRAPSCALE" in lnlst[0]: + lines.append(f"{lnlst[0]}") + else: + continue + + if openfile: + f.close() + + # reset unit numbers + unitnumber = cls._defaultunit() + if ext_unit_dict is not None: + for key, value in ext_unit_dict.items(): + if value.filetype == cls._ftype(): + unitnumber[0] = key + fname = os.path.basename(value.filename) + else: + fname = os.path.basename(filename) + + # initialize filenames list + filenames = [fname, None, None, None, None, None] + + # fill remainder of filenames list + if ihedun > 0: + unitnumber[1] = ihedun + try: + filenames[1] = os.path.basename(ext_unit_dict[ihedun].filename) + except: + if model.verbose: + print("head file name will be generated by flopy") + if iddnun > 0: + unitnumber[2] = iddnun + try: + filenames[2] = os.path.basename(ext_unit_dict[iddnun].filename) + except: + if model.verbose: + print("drawdown file name will be generated by flopy") + if ibouun > 0: + unitnumber[4] = ibouun + try: + filenames[4] = os.path.basename(ext_unit_dict[ibouun].filename) + except: + if model.verbose: + print("ibound file name will be generated by flopy") + if cboufm is None: + cboufm = True + if ispcun > 0: + unitnumber[5] = ispcun + try: + filenames[5] = os.path.basename(ext_unit_dict[ispcun].filename) + except: + if model.verbose: + print("concentration file name will be generated by flopy") + + # add unit numbers to pop_key_list + for u in unitnumber: + model.add_pop_key_list(u) + + return cls( + model, + ihedfm=ihedfm, + iddnfm=iddnfm, + ispcfm=ispcfm, + chedfm=chedfm, + cddnfm=cddnfm, + cboufm=cboufm, + cspcfm=cspcfm, + compact=compact, + stress_period_data=stress_period_data, + unitnumber=unitnumber, + filenames=filenames, + **kwargs, + ) + + @staticmethod + def _ftype(): + return "OC" + + @staticmethod + def _defaultunit(): + return [14, 0, 0, 0, 0, 0, 0] diff --git a/flopy/mfusg/mfusgpcb.py b/flopy/mfusg/mfusgpcb.py new file mode 100644 index 0000000000..19b1320c51 --- /dev/null +++ b/flopy/mfusg/mfusgpcb.py @@ -0,0 +1,294 @@ +""" +mfusgpcb module. Contains the MfUsgPcb class. Note that the user can access +the MfUsgPcb class as `flopy.mfusg.MfUsgPcb`. + +""" + +import numpy as np + +from ..pakbase import Package +from ..utils import MfList +from ..utils.recarray_utils import create_empty_recarray +from .mfusg import MfUsg + + +class MfUsgPcb(Package): + """ + MODFLOW USG Transport - Prescribed Concentration Boundary (PCB) Package Class. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to which + this package will be added. + ipakcb : int, optional + Toggles whether cell-by-cell budget data should be saved. If None or zero, + budget data will not be saved (default is None). + stress_period_data : list, recarray, dataframe or dictionary of boundaries. + Each pcb cell is defined through definition of + layer(int), row(int), column(int), stage(float), conductance(float) + The simplest form is a dictionary with a lists of boundaries for each + stress period, where each list of boundaries itself is a list of + boundaries. Indices of the dictionary are the numbers of the stress + period. This gives the form of:: + + stress_period_data = + {0: [ + [lay, row, col, iSpec, conc], + [lay, row, col, iSpec, conc], + [lay, row, col, iSpec, conc], + ], + 1: [ + [lay, row, col, iSpec, conc], + [lay, row, col, iSpec, conc], + [lay, row, col, iSpec, conc], + ], ... + kper: + [ + [lay, row, col, iSpec, conc], + [lay, row, col, iSpec, conc], + [lay, row, col, iSpec, conc], + ] + } + + Note that if no values are specified for a certain stress period, then + the list of boundaries for the previous stress period for which values + were defined is used. Full details of all options to specify + stress_period_data can be found in the flopy3boundaries Notebook in + the basic subdirectory of the examples directory + dtype : dtype definition + if data type is different from default + options : list of strings + Package options. (default is None). + extension : string + Filename extension (default is 'pcb') + unitnumber : int + File unit number (default is None). + filenames : str or list of str + Filenames to use for the package and the output files. If + filenames=None the package name will be created using the model name + and package extension and the cbc output name will be created using + the model name and .cbc extension (for example, modflowtest.cbc), + if ipakcb is a number greater than zero. If a single string is passed + the package will be set to the string and cbc output names will be + created using the model name and .cbc extension, if ipakcb is a + number greater than zero. To define the names for all package files + (input and output) the length of the list of strings should be 2. + Default is None. + + Attributes + ---------- + + Methods + ------- + + See Also + -------- + + Notes + ----- + Parameters are not supported in FloPy. + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg() + >>> lrcsc = {0:[2, 3, 4, 1, 100.]} #this pcb will be applied to all + >>> #stress periods + >>> pcb = flopy.mfusg.MfUsgPcb(ml, stress_period_data=lrcsc) + + """ + + def __init__( + self, + model, + ipakcb=None, + stress_period_data=None, + dtype=None, + no_print=False, + options=None, + extension="pcb", + unitnumber=None, + filenames=None, + ): + # set default unit number of one is not specified + if unitnumber is None: + unitnumber = MfUsgPcb._defaultunit() + + # set filenames + filenames = self._prepare_filenames(filenames, 2) + + # cbc output file + self.set_cbc_output_file(ipakcb, model, filenames[1]) + + # call base package constructor + super().__init__( + model, + extension=extension, + name=self._ftype(), + unit_number=unitnumber, + filenames=filenames[0], + ) + + self._generate_heading() + + self.no_print = no_print + self.np = 0 + if options is None: + options = [] + if self.no_print: + options.append("NOPRINT") + self.options = options + self.parent.add_package(self) + if dtype is not None: + self.dtype = dtype + else: + self.dtype = self.get_default_dtype(structured=self.parent.structured) + self.stress_period_data = MfList(self, stress_period_data) + + def _ncells(self): + """Maximum number of cells that have general head boundaries + (developed for MT3DMS SSM package). + + Returns + ------- + ncells: int + maximum number of pcb cells + + """ + return self.stress_period_data.mxact + + def write_file(self, check=False): + """ + Write the package file. + + Parameters + ---------- + check : boolean + Check package data for common errors. (default True) + + Returns + ------- + None + + """ + if check: # allows turning off package checks when writing files at model level + self.check( + f=f"{self.name[0]}.chk", + verbose=self.parent.verbose, + level=1, + ) + f_pcb = open(self.fn_path, "w") + f_pcb.write(f"{self.heading}\n") + f_pcb.write(f"{self.stress_period_data.mxact:10d}{self.ipakcb:10d}") + for option in self.options: + f_pcb.write(f" {option}") + f_pcb.write("\n") + self.stress_period_data.write_transient(f_pcb) + f_pcb.close() + + def add_record(self, kper, index, values): + try: + self.stress_period_data.add_record(kper, index, values) + except Exception as e: + raise Exception(f"MfUsgPcb error adding record to list: {e!s}") + + @staticmethod + def get_empty(ncells=0, aux_names=None, structured=True): + # get an empty recarray that corresponds to dtype + dtype = MfUsgPcb.get_default_dtype(structured=structured) + if aux_names is not None: + dtype = Package.add_to_dtype(dtype, aux_names, np.float32) + return create_empty_recarray(ncells, dtype, default_value=-1.0e10) + + @staticmethod + def get_default_dtype(structured=True): + if structured: + dtype = np.dtype( + [ + ("k", int), + ("i", int), + ("j", int), + ("iSpec", int), + ("conc", np.float32), + ] + ) + else: + dtype = np.dtype([("node", int), ("iSpec", int), ("conc", np.float32)]) + return dtype + + @staticmethod + def _get_sfac_columns(): + return ["conc"] + + @classmethod + def load(cls, f, model, nper=None, ext_unit_dict=None, check=False): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + nper : int + The number of stress periods. If nper is None, then nper will be + obtained from the model object. (default is None). + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + check : boolean + Check package data for common errors. (default True) + + Returns + ------- + pcb : MfUsgPcb object + + Examples + -------- + + >>> import flopy + >>> ml = flopy.mfusg.MfUsg(e) + >>> dis = flopy.modflow.ModflowDis.load('Test1.dis', ml) + >>> pcb = flopy.mfusg.MfUsgPcb.load('Test1.pcb', ml) + + """ + + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + if model.verbose: + print("loading pcb package file...") + + if model.version != "mfusg": + print( + "Warning: model version was reset from '{}' to 'mfusg' " + "in order to load a DDF file".format(model.version) + ) + model.version = "mfusg" + + return Package.load( + f, + model, + cls, + nper=nper, + check=check, + ext_unit_dict=ext_unit_dict, + ) + + @staticmethod + def _ftype(): + return "PCB" + + @staticmethod + def _defaultunit(): + return 154 diff --git a/flopy/mfusg/mfusgrch.py b/flopy/mfusg/mfusgrch.py new file mode 100644 index 0000000000..40bdedf937 --- /dev/null +++ b/flopy/mfusg/mfusgrch.py @@ -0,0 +1,518 @@ +""" +mfusgrch module. Contains the MfUsgRch class. Note that the user can access +the MfUsgRch class as `flopy.mfusg.MfUsgRch`. + +""" + +import numpy as np + +from ..modflow.mfparbc import ModflowParBc as mfparbc +from ..modflow.mfrch import ModflowRch +from ..utils import Transient2d, Transient3d, Util2d, read1d +from ..utils.flopy_io import line_parse +from ..utils.utils_def import ( + get_pak_vals_shape, + get_unitnumber_from_ext_unit_dict, +) +from .mfusg import MfUsg + + +class MfUsgRch(ModflowRch): + """ + MFUSG TRANSPORT Recharge Package Class. + + Parameters + ---------- + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to which + this package will be added. + ipakcb : int, optional + Toggles whether cell-by-cell budget data should be saved. If None or zero, + budget data will not be saved (default is None). + nrchop : int + is the recharge option code. + 1: Recharge to top grid layer only + 2: Recharge to layer defined in irch + 3: Recharge to highest active cell (default is 3). + rech : float or filename or ndarray or dict keyed on kper (zero-based) + Recharge flux (default is 1.e-3, which is used for all stress periods) + irch : int or filename or ndarray or dict keyed on kper (zero-based) + Layer (for an unstructured grid) or node (for an unstructured grid) to + which recharge is applied in each vertical column (only used when + nrchop=2). Default is 0, which is used for all stress periods. + extension : string + Filename extension (default is 'rch') + unitnumber : int + File unit number (default is None). + filenames : str or list of str + Filenames to use for the package and the output files. If + filenames=None the package name will be created using the model name + and package extension and the cbc output name will be created using + the model name and .cbc extension (for example, modflowtest.cbc), + if ipakcb is a number greater than zero. If a single string is passed + the package will be set to the string and cbc output names will be + created using the model name and .cbc extension, if ipakcb is a + number greater than zero. To define the names for all package files + (input and output) the length of the list of strings should be 2. + Default is None. + + Attributes + ---------- + + Methods + ------- + + See Also + -------- + + Notes + ----- + Parameters are supported in Flopy only when reading in existing models. + Parameter values are converted to native values in Flopy and the + connection to "parameters" is thus nonexistent. + + Examples + -------- + + """ + + def __init__( + self, + model, + nrchop=3, + ipakcb=None, + rech=1e-3, + irch=0, + seepelev=0, + mxrtzones=0, + iconc=0, + selev=None, + iznrch=None, + rchconc=None, + unitnumber=None, + filenames=None, + ): + """Package constructor.""" + msg = ( + "Model object must be of type flopy.mfusg.MfUsg\n" + f"but received type: {type(model)}." + ) + assert isinstance(model, MfUsg), msg + + # call base package constructor + super().__init__( + model, + nrchop=nrchop, + ipakcb=ipakcb, + rech=rech, + irch=irch, + unitnumber=unitnumber, + filenames=filenames, + ) + + self.mxrtzones = mxrtzones + self.seepelev = seepelev + self.iconc = iconc + + self.selev = None + if selev is not None: + selev_u2d_shape = get_pak_vals_shape(model, selev) + self.selev = Transient2d( + model, selev_u2d_shape, np.float32, rech, name="rech_selev" + ) + + self.iznrch = None + if iznrch is not None: + iznrch_u2d_shape = get_pak_vals_shape(model, iznrch) + self.iznrch = Transient2d( + model, iznrch_u2d_shape, np.int32, rech, name="rech_izn" + ) + + self.rchconc = None + if rchconc is not None: + mcomp = model.mcomp + if model.iheat > 0: + mcomp = model.mcomp + 1 + self.irchconc = [1] * mcomp + rchconc_u3d_shape = (mcomp,) + get_pak_vals_shape(model, rech) + self.rchconc = Transient3d( + model, rchconc_u3d_shape, np.float32, rchconc, name="rech_conc" + ) + + def write_file(self, check=True, f=None): + """ + Write the package file. + + Parameters + ---------- + check : boolean + Check package data for common errors. (default True) + + Returns + ------- + None + + """ + # allows turning off package checks when writing files at model level + if check: + self.check( + f=f"{self.name[0]}.chk", + verbose=self.parent.verbose, + level=1, + ) + nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper + # Open file for writing + if f is not None: + f_rch = f + else: + f_rch = open(self.fn_path, "w") + f_rch.write(f"{self.heading}\n") + f_rch.write(f"{self.nrchop:10d}{self.ipakcb:10d}") + if self.seepelev: + f_rch.write(" SEEPELEV") + if self.iconc: + f_rch.write(" CONC") + if self.mxrtzones: + f_rch.write(f" RTS {self.mxrtzones:4.0d}") + f_rch.write("\n") + + mcomp = self.parent.mcomp + if self.parent.iheat > 0: + mcomp = mcomp + 1 + + if self.iconc: + for icomp in range(mcomp): + f_rch.write(f"{self.irchconc[icomp]:10.0f}") + f_rch.write("\n") + + if self.nrchop == 2: + irch = {} + for kper, u2d in self.irch.transient_2ds.items(): + irch[kper] = u2d.array + 1 + irch = Transient2d( + self.parent, + self.irch.shape, + self.irch.dtype, + irch, + self.irch.name, + ) + if not self.parent.structured: + mxndrch = np.max( + [u2d.array.size for kper, u2d in self.irch.transient_2ds.items()] + ) + f_rch.write(f"{mxndrch:10d}\n") + + for kper in range(nper): + inrech, file_entry_rech = self.rech.get_kper_entry(kper) + if self.nrchop == 2: + inirch, file_entry_irch = irch.get_kper_entry(kper) + if not self.parent.structured: + inirch = self.rech[kper].array.size + else: + inirch = -1 + + f_rch.write(f"{inrech:10d}{inirch:10d} ") + + if self.iznrch is not None: + f_rch.write(" INRCHZONES 1") + if self.selev is not None: + f_rch.write(" INSELEV 1") + if self.rchconc is not None: + f_rch.write(" INCONC 1") + f_rch.write("# Stress period {kper + 1}\n") + + if inrech >= 0: + f_rch.write(file_entry_rech) + if self.nrchop == 2: + if inirch >= 0: + f_rch.write(file_entry_irch) + + if self.iznrch is not None: + inrech, file_entry_iznrch = self.iznrch.get_kper_entry(kper) + f_rch.write(file_entry_iznrch) + if self.selev is not None: + inrech, file_entry_selev = self.selev.get_kper_entry(kper) + f_rch.write(file_entry_selev) + if self.rchconc is not None: + inrech, file_entry_rchconc = self.rchconc.get_kper_entry(kper) + f_rch.write(file_entry_rchconc) + f_rch.close() + + @classmethod + def load(cls, f, model, nper=None, ext_unit_dict=None, check=True): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + nper : int + The number of stress periods. If nper is None, then nper will be + obtained from the model object. (default is None). + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + check : boolean + Check package data for common errors. (default True) + + Returns + ------- + rch : MfUsgRch object + MfUsgRch object. + + Examples + -------- + + >>> import flopy + >>> m = flopy.modflow.Modflow() + >>> rch = flopy.modflow.MfUsgRch.load('test.rch', m) + + """ + if model.verbose: + print("loading rch package file...") + + openfile = not hasattr(f, "read") + if openfile: + filename = f + f = open(filename, "r") + + # dataset 0 -- header + while True: + line = f.readline() + if line[0] != "#": + break + npar = 0 + + if "parameter" in line.lower(): + raw = line.strip().split() + npar = int(raw[1]) + if npar > 0: + if model.verbose: + print(f" Parameters detected. Number of parameters = {npar}") + line = f.readline() + # dataset 2 + t = line_parse(line) + nrchop = int(t[0]) + ipakcb = int(t[1]) + + # item 2 - options + seepelev = 0 + if "SEEPELEV" in t: + seepelev = 1 + + iconc = 0 + if "CONC" in line: + iconc = 1 + + mxrtzones = 0 + if "RTS" in t: + idx = t.index("RTS") + mxrtzones = int(t[idx + 1]) + + # dataset 2b for mfusg + if not model.structured and nrchop == 2: + line = f.readline() + t = line_parse(line) + mxndrch = int(t[0]) + + mcomp = model.mcomp + if model.iheat > 0: + mcomp = model.mcomp + 1 + + irchconc = np.empty((mcomp), dtype=np.int32) + if iconc: + if model.verbose: + print(f" loading IRCHCONC[{mcomp}] ...") + irchconc = read1d(f, irchconc) + + print(f"irchconc: {irchconc}") + + # dataset 3 and 4 - parameters data + pak_parms = None + if npar > 0: + pak_parms = mfparbc.loadarray(f, npar, model.verbose) + + if nper is None: + nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper() + else: + nrow, ncol, nlay, _ = model.get_nrow_ncol_nlay_nper() + + u2d_shape = (nrow, ncol) + + # read data for every stress period + current_rech = [] + rech = {} + irch = None + if nrchop == 2: + irch = {} + current_irch = [] + + iznrch = None + if mxrtzones: + iznrch = {} + iniznrch = [0] * nper + current_iznrch = [] + + selev = None + if seepelev: + selev = {} + inselev = [0] * nper + current_selev = [] + + rchconc = None + if iconc: + rchconc = {} + inconc = [0] * nper + current_rchconc = [] + + for iper in range(nper): + line = f.readline() + t = line_parse(line) + inrech = int(t[0]) + + if nrchop == 2: + inirch = int(t[1]) + if (not model.structured) and (inirch >= 0): + u2d_shape = (1, inirch) + elif not model.structured: + u2d_shape = (1, ncol[0]) + + if "INSELEV" in t: + idx = t.index("INSELEV") + inselev[iper] = int(t[idx + 1]) + + if "INRCHZONES" in t: + idx = t.index("INIZNRCH") + iniznrch[iper] = int(t[idx + 1]) + + if "INCONC" in t: + inconc[iper] = 1 + + if inrech >= 0: + if npar == 0: + if model.verbose: + print(f" loading rech stress period {iper + 1:3d}...") + t = Util2d.load( + f, + model, + u2d_shape, + np.float32, + "rech", + ext_unit_dict, + ) + else: + parm_dict = {} + for ipar in range(inrech): + line = f.readline() + t = line.strip().split() + pname = t[0].lower() + try: + c = t[1].lower() + instance_dict = pak_parms.bc_parms[pname][1] + if c in instance_dict: + iname = c + else: + iname = "static" + except: + iname = "static" + parm_dict[pname] = iname + t = mfparbc.parameter_bcfill(model, u2d_shape, parm_dict, pak_parms) + + current_rech = t + rech[iper] = current_rech + + if nrchop == 2: + if inirch >= 0: + if model.verbose: + print(f" loading irch stress period {iper + 1:3d}...") + t = Util2d.load( + f, model, u2d_shape, np.int32, "irch", ext_unit_dict + ) + current_irch = Util2d( + model, u2d_shape, np.int32, t.array - 1, "irch" + ) + irch[iper] = current_irch + + # Item 9 for mfusg transport + if mxrtzones: + if iniznrch[iper]: + if model.verbose: + print(f" loading iznrch stress period {iper + 1:3d}...") + current_iznrch = Util2d.load( + f, model, u2d_shape, np.int32, "iznrch", ext_unit_dict + ) + iznrch[iper] = current_iznrch + + # Item 10 for mfusg transport + if seepelev: + if inselev[iper]: + if model.verbose: + print(f" loading selev stress period {iper + 1:3d}...") + current_selev = Util2d.load( + f, model, u2d_shape, np.float32, "selev", ext_unit_dict + ) + selev[iper] = current_selev + + # Item 11 for mfusg transport + if iconc: + if inconc[iper]: + if model.verbose: + print(f" loading rch conc stress period {iper + 1:3d}...") + if model.iheat > 0: + mcomp = model.mcomp + 1 + + current_rchconc = [0] * mcomp + for icomp in range(mcomp): + if irchconc[icomp]: + if model.verbose: + print( + f" loading rch conc stress period {iper + 1:3d}" + f" component {icomp + 1}..." + ) + current_rchconc[icomp] = Util2d.load( + f, + model, + u2d_shape, + np.float32, + "rchconc", + ext_unit_dict, + ) + rchconc[iper] = current_rchconc + + if openfile: + f.close() + + unitnumber, filenames = get_unitnumber_from_ext_unit_dict( + model, cls, ext_unit_dict, ipakcb + ) + + # create recharge package instance + rch = cls( + model, + nrchop=nrchop, + ipakcb=ipakcb, + rech=rech, + irch=irch, + seepelev=seepelev, + mxrtzones=mxrtzones, + iconc=iconc, + selev=selev, + iznrch=iznrch, + rchconc=rchconc, + unitnumber=unitnumber, + filenames=filenames, + ) + if check: + rch.check( + f=f"{rch.name[0]}.chk", + verbose=rch.parent.verbose, + level=0, + ) + return rch diff --git a/flopy/mfusg/mfusgsms.py b/flopy/mfusg/mfusgsms.py index afd00d11c2..c8fb2619eb 100644 --- a/flopy/mfusg/mfusgsms.py +++ b/flopy/mfusg/mfusgsms.py @@ -313,12 +313,21 @@ def __init__( self.rrctol = rrctol self.idroptol = idroptol self.epsrn = epsrn - self.clin = clin + + if clin is None: + self.clin = "CG" + else: + self.clin = clin + + if relaxpcgu is None: + self.relaxpcgu = 0.97 + else: + self.relaxpcgu = relaxpcgu + self.ipc = ipc self.iscl = iscl self.iord = iord self.rclosepcgu = rclosepcgu - self.relaxpcgu = relaxpcgu if options is None: self.options = [] else: @@ -343,52 +352,27 @@ def write_file(self): if nopt > 0: f.write(" ".join(self.options) + "\n") f.write( - "{} {} {} {} {} {} {}\n".format( - self.hclose, - self.hiclose, - self.mxiter, - self.iter1, - self.iprsms, - self.nonlinmeth, - self.linmeth, - ) + f" {self.hclose:9.2e} {self.hiclose:9.2e} {self.mxiter:9d} {self.iter1:9d}" + f" {self.iprsms:9d} {self.nonlinmeth:9d} {self.linmeth:9d}\n" ) + if self.nonlinmeth != 0 and nopt == 0: f.write( - "{} {} {} {} {} {} {} {}\n".format( - self.theta, - self.akappa, - self.gamma, - self.amomentum, - self.numtrack, - self.btol, - self.breduc, - self.reslim, - ) + f" {self.theta:9.2e} {self.akappa:9.2e} {self.gamma:9.2e}" + f" {self.amomentum:9.2e} {self.numtrack:9d} {self.btol:9.2e}" + f" {self.breduc:9.2e} {self.reslim:9.2e}\n" ) + if self.linmeth == 1 and nopt == 0: f.write( - "{} {} {} {} {} {} {} {}\n".format( - self.iacl, - self.norder, - self.level, - self.north, - self.iredsys, - self.rrctol, - self.idroptol, - self.epsrn, - ) + f" {self.iacl:9d} {self.norder:9d} {self.level:9d} {self.north:9d}" + f" {self.iredsys:9d} {self.rrctol:9.2e} {self.idroptol:9d}" + f" {self.epsrn:9.2e}\n" ) if self.linmeth == 2 and nopt == 0: f.write( - "{} {} {} {} {} {}\n".format( - self.clin, - self.ipc, - self.iscl, - self.iord, - self.rclosepcgu, - self.relaxpcgu, - ) + f" {self.clin:9s} {self.ipc:9d} {self.iscl:9d} {self.iord:9d}" + f" {self.rclosepcgu:9.2e} {self.relaxpcgu:9.2e}\n" ) f.write("\n") f.close() @@ -473,13 +457,10 @@ def load(cls, f, model, ext_unit_dict=None): nonlinmeth = int(ll.pop(0)) linmeth = int(ll.pop(0)) if model.verbose: - print(f" HCLOSE {hclose}") - print(f" HICLOSE {hiclose}") - print(f" MXITER {mxiter}") - print(f" ITER1 {iter1}") - print(f" IPRSMS {iprsms}") - print(f" NONLINMETH {nonlinmeth}") - print(f" LINMETH {linmeth}") + print( + f" HCLOSE {hclose} HICLOSE {hiclose} MXITER {mxiter} ITER1 {iter1}" + f" IPRSMS {iprsms} NONLINMETH {nonlinmeth} LINMETH {linmeth}" + ) # Record 2 theta = None @@ -510,14 +491,10 @@ def load(cls, f, model, ext_unit_dict=None): breduc = float(ll.pop(0)) reslim = float(ll.pop(0)) if model.verbose: - print(f" THETA {theta}") - print(f" AKAPPA {akappa}") - print(f" GAMMA {gamma}") - print(f" AMOMENTUM {amomentum}") - print(f" NUMTRACK {numtrack}") - print(f" BTOL {btol}") - print(f" BREDUC {breduc}") - print(f" RESLIM {reslim}") + print( + f" THETA {theta} AKAPPA{akappa} GAMMA{gamma} AMOMENTUM{amomentum}" + f" NUMTRACK {numtrack} BTOL {btol} BREDUC {breduc} RESLIM {reslim}" + ) iacl = None norder = None @@ -547,14 +524,11 @@ def load(cls, f, model, ext_unit_dict=None): idroptol = int(ll.pop(0)) epsrn = float(ll.pop(0)) if model.verbose: - print(f" IACL {iacl}") - print(f" NORDER {norder}") - print(f" LEVEL {level}") - print(f" NORTH {north}") - print(f" IREDSYS {iredsys}") - print(f" RRCTOL {rrctol}") - print(f" IDROPTOL {idroptol}") - print(f" EPSRN {epsrn}") + print( + f" IACL {iacl} NORDER {norder} LEVEL {level} NORTH {north}" + f" IREDSYS{iredsys} RRCTOL{rrctol} IDROPTOL{idroptol}" + f" EPSRN{epsrn}" + ) clin = None ipc = None @@ -579,12 +553,10 @@ def load(cls, f, model, ext_unit_dict=None): if len(ll) > 0: relaxpcgu = float(ll.pop(0)) if model.verbose: - print(f" CLIN {clin}") - print(f" IPC {ipc}") - print(f" ISCL {iscl}") - print(f" IORD {iord}") - print(f" RCLOSEPCGU {rclosepcgu}") - print(f" RELAXPCGU {relaxpcgu}") + print( + f" CLIN {clin} IPC {ipc} ISCL {iscl} IORD {iord} " + f" RCLOSEPCGU {rclosepcgu} RELAXPCGU {relaxpcgu}" + ) if openfile: f.close() diff --git a/flopy/mfusg/mfusgwel.py b/flopy/mfusg/mfusgwel.py index 9a859734b4..8fdb47fb56 100644 --- a/flopy/mfusg/mfusgwel.py +++ b/flopy/mfusg/mfusgwel.py @@ -3,10 +3,6 @@ Contains the MfUsgWel class. Note that the user can access the MfUsgWel class as `flopy.mfusg.MfUsgWel`. - -Additional information for this MODFLOW package can be found at the `Online -MODFLOW Guide -`_. """ from ..modflow.mfwel import ModflowWel diff --git a/flopy/modflow/mffhb.py b/flopy/modflow/mffhb.py index 259caa8ba7..0c29bf5482 100644 --- a/flopy/modflow/mffhb.py +++ b/flopy/modflow/mffhb.py @@ -623,7 +623,10 @@ def load(cls, f, model, nper=None, ext_unit_dict=None): structured=model.structured, ) for n in range(nhed): - tds7 = read1d(f, np.empty((nbdtim + 4,))) + if model.structured: + tds7 = read1d(f, np.empty((nbdtim + 4,))) + else: + tds7 = read1d(f, np.empty((nbdtim + 2,))) ds7[n] = tuple(tds7) if model.structured: