"""
mfdisu module. Contains the MfUsgDisU class. Note that the user can access
the MfUsgDisU class as `flopy.mfusg.MfUsgDisU`.
"""
import numpy as np
from ..discretization.unstructuredgrid import UnstructuredGrid
from ..pakbase import Package
from ..utils import Util2d, Util3d, read1d
from ..utils.reference import TemporalReference
from .mfusg import MfUsg
ITMUNI = {"u": 0, "s": 1, "m": 2, "h": 3, "d": 4, "y": 5}
LENUNI = {"u": 0, "f": 1, "m": 2, "c": 3}
[docs]class MfUsgDisU(Package):
"""
MODFLOW Unstructured Discretization Package Class.
Parameters
----------
model : model object
The model object (of type :class:`flopy.mfusg.MfUsg`) to which
this package will be added.
nodes : int
Number of nodes in the model grid (default is 2).
nlay : int
Number of layers in the model grid (default is 1).
njag : int
Total number of connections of an unstructured grid. njag is used to
dimension the sparse matrix in a compressed row storage format. For
symmetric arrays, only the upper triangle of the matrix may be
entered. For that case, the symmetric portion (minus the diagonal
terms) is dimensioned as njags = (njag - nodes) / 2.
(default is None).
ivsd : int
is the vertical sub-discretization index. For purposes of this flag,
vertical sub-discretization is defined to occur when all layers are
not a stacked representation of each other.
If IVSD = 0 there is no sub-discretization of layers within the model
domain. That is, grids are not nested in the vertical direction.
However, one layer may have a different grid structure from the next
due to different sub-gridding structures within each layer.
If IVSD = 1 there could be sub-discretization of layers with
vertically nested grids (as shown in Figure 5c in the MODFLOW-USG
document) within the domain. For this case, the vertical connection
index IVC is required to determine the vertical connections of every
node. Otherwise, the vertical connections are internally computed and
IVC is not read.
If IVSD = -1 there is no vertical sub-discretization of layers, and
further, the horizontal discretization of all layers is the same. For
this case, the cell areas (AREA) are read only for one layer and are
computed to be the same for all the stacked layers. A structured
finite-difference grid is an example of this condition.
(default is 0).
nper : int
Number of model stress periods (the default is 1).
itmuni : int
Time units, default is days (4)
lenuni : int
Length units, default is meters (2)
idsymrd : int
A flag indicating if the finite-volume connectivity information of an
unstructured grid is input as a full matrix or as a symmetric matrix
in the input file.
If idsymrd is 0 the finite-volume connectivity information is provided
for the full matrix of the porous matrix grid-block connections of an
unstructured grid. The code internally stores only the symmetric
portion of this information. This input structure (IDSYMRD=0) is easy
to organize but contains unwanted information which is parsed out
when the information is stored.
If idsymrd is 1 then finite-volume connectivity information is
provided only for the upper triangular portion of the porous matrix
grid-block connections within the unstructured grid. This input
structure (IDSYMRD=1) is compact but is slightly more complicated to
organize. Only the non-zero upper triangular items of each row are
read in sequence for all symmetric matrices.
(default is 0).
laycbd : int or array of ints (nlay), optional
An array of flags indicating whether or not a layer has a Quasi-3D
confining bed below it. 0 indicates no confining bed, and not zero
indicates a confining bed. LAYCBD for the bottom layer must be 0. (the
default is 0)
nodelay : int or array of ints (nlay)
The number of cells in each layer. (the default is None, which means
the number of cells in a layer is equal to nodes / nlay).
top : float or array of floats (nodes), optional
An array of the top elevation for every cell. For the situation in
which the top layer represents a water-table aquifer, it may be
reasonable to set Top equal to land-surface elevation (the default is
1.0)
bot : float or array of floats (nodes), optional
An array of the bottom elevation for each model cell (the default is
0.)
area : float or array of floats
Surface area for model cells. Area is for only one layer if IVSD = -1
to indicate that the grid is vertically stacked. Otherwise, area is
required for each layer in the model grid. Note that there may be
different number of nodes per layer (ndslay) for an unstructured grid.
(default is 1.0)
iac : array of integers
is a vector indicating the number of connections plus 1 for each
node. Note that the IAC array is only supplied for the GWF cells;
the IAC array is internally expanded to include CLN or GNC nodes if
they are present in a simulation.
(default is None. iac must be provided).
ja : array of integers
is a list of cell number (n) followed by its connecting cell numbers
(m) for each of the m cells connected to cell n. This list is
sequentially provided for the first to the last GWF cell. Note that
the cell and its connections are only supplied for the GWF cells and
their connections to the other GWF cells. This connectivity is
internally expanded if CLN or GNC nodes are present in a simulation.
Also note that the JA list input may be chopped up to have every node
number and its connectivity list on a separate line for ease in
readability of the file. To further ease readability of the file, the
node number of the cell whose connectivity is subsequently listed,
may be expressed as a negative number the sign of which is
subsequently corrected by the code.
(default is None. ja must be provided).
ivc : int or array of integers
is an index array indicating the direction between a node n and all
its m connections. IVC = 0 if the connection between n and m is
horizontal. IVC = 1 if the connecting node m is vertically oriented
to node n. Note that if the CLN Process is active, the connection
between two CLN cells has IVC = 2 and the connection between a CLN
cell and a GWF cell has IVC = 3.
(default is None. ivc must be provided if ivsd = 1)
cl1 : float or array of floats
is the perpendicular length between the center of a node (node 1) and
the interface between the node and its adjoining node (node 2).
(default is None. cl1 and cl2 must be specified, or cl12 must be
specified)
cl2 : float or array of floats
is the perpendicular length between node 2 and the interface between
nodes 1 and 2, and is at the symmetric location of CL1.
(default is None. cl1 and cl2 must be specified, or cl12 must be
specified)
cl12 : float or array of floats
is the array containing CL1 and CL2 lengths, where CL1 is the
perpendicular length between the center of a node (node 1) and the
interface between the node and its adjoining node (node 2). CL2,
which is the perpendicular length between node 2 and the interface
between nodes 1 and 2 is at the symmetric location of CL1. The array
CL12 reads both CL1 and CL2 in the upper and lower triangular
portions of the matrix respectively. Note that the CL1 and CL2 arrays
are only supplied for the GWF cell connections and are internally
expanded if CLN or GNC nodes exist in a simulation.
(default is None. cl1 and cl2 must be specified, or cl12 must be
specified)
fahl : float or array of floats
Area of the interface Anm between nodes n and m.
(default is None. fahl must be specified.)
perlen : float or array of floats (nper)
An array of the stress period lengths.
nstp : int or array of ints (nper)
Number of time steps in each stress period (default is 1).
tsmult : float or array of floats (nper)
Time step multiplier (default is 1.0).
steady : bool or array of bool (nper)
True or False indicating whether or not stress period is steady state
(default is True).
extension : string
Filename extension (default is 'dis')
unitnumber : int
File unit number (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
----------
heading : str
Text string written to top of package input file.
Methods
-------
See Also
--------
Notes
-----
Now works for multi-layer USG models since u3d was modified to handle
multiple u2d instances of different size.
Examples
--------
>>> import flopy
>>> m = flopy.mfusg.MfUsg()
>>> disu = flopy.mfusg.MfUsgDisU(m)
"""
def __init__(
self,
model,
nodes=2,
nlay=1,
njag=None,
ivsd=0,
nper=1,
itmuni=4,
lenuni=2,
idsymrd=0,
laycbd=0,
nodelay=None,
top=1,
bot=0,
area=1.0,
iac=None,
ja=None,
ivc=None,
cl1=None,
cl2=None,
cl12=None,
fahl=None,
perlen=1,
nstp=1,
tsmult=1,
steady=True,
extension="disu",
unitnumber=None,
filenames=None,
start_datetime=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),
)
# Set values of all parameters
self._generate_heading()
self.nodes = nodes
self.nlay = nlay
self.njag = njag
self.ivsd = ivsd
self.nper = nper
try:
self.itmuni = int(itmuni)
except:
self.itmuni = ITMUNI[itmuni.lower()[0]]
try:
self.lenuni = int(lenuni)
except:
self.lenuni = LENUNI[lenuni.lower()[0]]
self.idsymrd = idsymrd
# LAYCBD
self.laycbd = Util2d(
model, (self.nlay,), np.int32, laycbd, name="laycbd"
)
self.laycbd[-1] = 0 # bottom layer must be zero
# NODELAY
if nodelay is None:
npl = int(nodes / nlay)
nodelay = []
for k in range(self.nlay):
nodelay.append(npl)
self.nodelay = Util2d(
model,
(self.nlay,),
np.int32,
nodelay,
name="nodelay",
locat=self.unit_number[0],
)
# set ncol and nrow for array readers
nrow = None
ncol = self.nodelay.array[:]
# Top and bot are both 1d arrays of size nodes
self.top = Util3d(
model,
(nlay, nrow, ncol),
np.float32,
top,
name="top",
locat=self.unit_number[0],
)
self.bot = Util3d(
model,
(nlay, nrow, ncol),
np.float32,
bot,
name="bot",
locat=self.unit_number[0],
)
# Area is Util2d if ivsd == -1, otherwise it is Util3d
if ivsd == -1:
self.area = Util2d(
model,
(self.nodelay[0],),
np.float32,
area,
"area",
locat=self.unit_number[0],
)
else:
self.area = Util3d(
model,
(nlay, nrow, ncol),
np.float32,
area,
name="area",
locat=self.unit_number[0],
)
# Connectivity and ivc
if iac is None:
raise Exception("iac must be provided")
self.iac = Util2d(
model,
(self.nodes,),
np.int32,
iac,
name="iac",
locat=self.unit_number[0],
)
assert self.iac.array.sum() == njag, "The sum of iac must equal njag."
if ja is None:
raise Exception("ja must be provided")
if ja[0] == 0:
# convert from zero-based to one-based
ja += 1
self.ja = Util2d(
model,
(self.njag,),
np.int32,
ja,
name="ja",
locat=self.unit_number[0],
)
self.ivc = None
if self.ivsd == 1:
if ivc is None:
raise Exception("ivc must be provided if ivsd is 1.")
self.ivc = Util2d(
model,
(self.njag,),
np.int32,
ivc,
name="ivc",
locat=self.unit_number[0],
)
# Connection lengths
if idsymrd == 1:
njags = int((njag - nodes) / 2)
if cl1 is None:
raise Exception("idsymrd is 1 but cl1 was not specified.")
if cl2 is None:
raise Exception("idsymrd is 1 but cl2 was not specified.")
self.cl1 = Util2d(
model,
(njags,),
np.float32,
cl1,
name="cl1",
locat=self.unit_number[0],
)
self.cl2 = Util2d(
model,
(njags,),
np.float32,
cl2,
name="cl2",
locat=self.unit_number[0],
)
if idsymrd == 0:
if cl12 is None:
raise Exception("idsymrd is 0 but cl12 was not specified")
self.cl12 = Util2d(
model,
(self.njag,),
np.float32,
cl12,
name="cl12",
locat=self.unit_number[0],
)
# Flow area (set size of array to njag or njags depending on idsymrd)
if fahl is None:
raise Exception("fahl must be provided")
if idsymrd == 1:
n = njags
elif idsymrd == 0:
n = self.njag
self.fahl = Util2d(
model,
(n,),
np.float32,
fahl,
name="fahl",
locat=self.unit_number[0],
)
# Stress period information
self.perlen = Util2d(
model, (self.nper,), np.float32, perlen, name="perlen"
)
self.nstp = Util2d(model, (self.nper,), np.int32, nstp, name="nstp")
self.tsmult = Util2d(
model, (self.nper,), np.float32, tsmult, name="tsmult"
)
self.steady = Util2d(model, (self.nper,), bool, steady, name="steady")
self.itmuni_dict = {
0: "undefined",
1: "seconds",
2: "minutes",
3: "hours",
4: "days",
5: "years",
}
if start_datetime is None:
start_datetime = model._start_datetime
if model.modelgrid is None:
model.modelgrid = UnstructuredGrid(
ncpl=self.nodelay.array,
top=self.top.array,
botm=self.bot.array,
lenuni=self.lenuni,
)
self.tr = TemporalReference(
itmuni=self.itmuni, start_datetime=start_datetime
)
self.start_datetime = start_datetime
# get neighboring nodes
self._get_neighboring_nodes()
# Add package and return
self.parent.add_package(self)
return
[docs] def checklayerthickness(self):
"""
Check layer thickness.
"""
return (self.parent.modelgrid.cell_thickness > 0).all()
[docs] def get_cell_volumes(self):
"""
Get an array of cell volumes.
Returns
-------
vol : array of floats (nodes)
"""
vol = np.empty((self.nodes,))
for n in range(self.nodes):
nn = n
if self.ivsd == -1:
nn = n % self.nodelay[0]
area = self.area[nn]
vol[n] = area * (self.top[n] - self.bot[n])
return vol
@property
def zcentroids(self):
"""
Return an array of size nodes that contains the vertical cell center
elevation.
"""
z = np.empty((self.nodes,))
z[:] = (self.top.array - self.bot.array) / 2.0
return z
@property
def ncpl(self):
return self.nodes / self.nlay
[docs] @classmethod
def load(cls, f, model, 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.mfusg.MfUsg`) 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`.
check : bool
Check package data for common errors. (default False; not setup yet)
Returns
-------
dis : MfUsgDisU object
MfUsgDisU object.
Examples
--------
>>> import flopy
>>> m = flopy.mfusg.MfUsg()
>>> disu = flopy.mfusg.MfUsgDisU.load('test.disu', m)
"""
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 disu package file...")
if model.version != "mfusg":
print(
"Warning: model version was reset from '{}' to 'mfusg' "
"in order to load a DISU file".format(model.version)
)
model.version = "mfusg"
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
# dataset 1
if model.verbose:
print(
" loading NODES, NLAY, NJAG, IVSD, NPER, ITMUNI, LENUNI,"
" IDSYMRD..."
)
ll = line.strip().split()
nodes = int(ll.pop(0))
nlay = int(ll.pop(0))
njag = int(ll.pop(0))
ivsd = int(ll.pop(0))
nper = int(ll.pop(0))
# mimic urword behavior in case these values aren't present on line
if len(ll) > 0:
itmuni = int(ll.pop(0))
else:
itmuni = 0
if len(ll) > 0:
lenuni = int(ll.pop(0))
else:
lenuni = 0
if len(ll) > 0:
idsymrd = int(ll.pop(0))
else:
idsymrd = 0
if model.verbose:
print(f" NODES {nodes}")
print(f" NLAY {nlay}")
print(f" NJAG {njag}")
print(f" IVSD {ivsd}")
print(f" NPER {nper}")
print(f" ITMUNI {itmuni}")
print(f" LENUNI {lenuni}")
print(f" IDSYMRD {idsymrd}")
# Calculate njags
njags = int((njag - nodes) / 2)
if model.verbose:
print(f" NJAGS calculated as {njags}")
# dataset 2 -- laycbd
if model.verbose:
print(" loading LAYCBD...")
laycbd = np.empty((nlay,), np.int32)
laycbd = read1d(f, laycbd)
if model.verbose:
print(f" LAYCBD {laycbd}")
# dataset 3 -- nodelay
if model.verbose:
print(" loading NODELAY...")
nodelay = Util2d.load(
f, model, (nlay,), np.int32, "nodelay", ext_unit_dict
)
if model.verbose:
print(f" NODELAY {nodelay}")
# dataset 4 -- top
if model.verbose:
print(" loading TOP...")
top = [0] * nlay
for k in range(nlay):
tpk = Util2d.load(
f, model, (nodelay[k],), np.float32, "top", ext_unit_dict
)
top[k] = tpk
if model.verbose:
for k, tpk in enumerate(top):
print(f" TOP layer {k}: {tpk.array}")
# dataset 5 -- bot
if model.verbose:
print(" loading BOT...")
bot = [0] * nlay
for k in range(nlay):
btk = Util2d.load(
f, model, (nodelay[k],), np.float32, "btk", ext_unit_dict
)
bot[k] = btk
if model.verbose:
for k, btk in enumerate(bot):
print(f" BOT layer {k}: {btk.array}")
# dataset 6 -- area
if model.verbose:
print(" loading AREA...")
if ivsd == -1:
area = Util2d.load(
f, model, (nodelay[0],), np.float32, "area", ext_unit_dict
)
else:
area = [0] * nlay
for k in range(nlay):
ak = Util2d.load(
f, model, (nodelay[k],), np.float32, "ak", ext_unit_dict
)
area[k] = ak
if model.verbose:
for k, ak in enumerate(area):
print(f" AREA layer {k}: {ak}")
# dataset 7 -- iac
if model.verbose:
print(" loading IAC...")
iac = Util2d.load(f, model, (nodes,), np.int32, "iac", ext_unit_dict)
if model.verbose:
print(f" IAC {iac}")
# dataset 8 -- ja
if model.verbose:
print(" loading JA...")
ja = Util2d.load(f, model, (njag,), np.int32, "ja", ext_unit_dict)
if model.verbose:
print(f" JA {ja}")
# dataset 9 -- ivc
ivc = None
if ivsd == 1:
if model.verbose:
print(" loading IVC...")
ivc = Util2d.load(
f, model, (njag,), np.int32, "ivc", ext_unit_dict
)
if model.verbose:
print(f" IVC {ivc}")
# dataset 10a -- cl1
cl1 = None
if idsymrd == 1:
if model.verbose:
print(" loading CL1...")
cl1 = Util2d.load(
f, model, (njags,), np.float32, "cl1", ext_unit_dict
)
if model.verbose:
print(f" CL1 {cl1}")
# dataset 10b -- cl2
cl2 = None
if idsymrd == 1:
if model.verbose:
print(" loading CL2...")
cl2 = Util2d.load(
f, model, (njags,), np.float32, "cl2", ext_unit_dict
)
if model.verbose:
print(f" CL2 {cl2}")
# dataset 11 -- cl12
cl12 = None
if idsymrd == 0:
if model.verbose:
print(" loading CL12...")
cl12 = Util2d.load(
f, model, (njag,), np.float32, "cl12", ext_unit_dict
)
if model.verbose:
print(f" CL12 {cl12}")
# dataset 12 -- fahl
fahl = None
if idsymrd == 0:
n = njag
elif idsymrd == 1:
n = njags
if model.verbose:
print(" loading FAHL...")
fahl = Util2d.load(f, model, (n,), np.float32, "fahl", ext_unit_dict)
if model.verbose:
print(f" FAHL {fahl}")
# dataset 7 -- stress period info
if model.verbose:
print(" loading stress period data...")
perlen = []
nstp = []
tsmult = []
steady = []
for k in range(nper):
line = f.readline()
a1, a2, a3, a4 = line.strip().split()[0:4]
a1 = float(a1)
a2 = int(a2)
a3 = float(a3)
if a4.upper() == "TR":
a4 = False
else:
a4 = True
perlen.append(a1)
nstp.append(a2)
tsmult.append(a3)
steady.append(a4)
if model.verbose:
print(f" PERLEN {perlen}")
print(f" NSTP {nstp}")
print(f" TSMULT {tsmult}")
print(f" STEADY {steady}")
if openfile:
f.close()
# set package unit number
unitnumber = None
filenames = [None]
if ext_unit_dict is not None:
unitnumber, filenames[0] = model.get_ext_dict_attr(
ext_unit_dict, filetype=cls._ftype()
)
# create dis object instance
disu = cls(
model,
nodes=nodes,
nlay=nlay,
njag=njag,
ivsd=ivsd,
nper=nper,
itmuni=itmuni,
lenuni=lenuni,
idsymrd=idsymrd,
laycbd=laycbd,
nodelay=nodelay,
top=top,
bot=bot,
area=area,
iac=iac,
ja=ja,
ivc=ivc,
cl1=cl1,
cl2=cl2,
cl12=cl12,
fahl=fahl,
perlen=perlen,
nstp=nstp,
tsmult=tsmult,
steady=steady,
unitnumber=unitnumber,
filenames=filenames,
)
# return dis object instance
return disu
[docs] def write_file(self):
"""
Write the package file.
Returns
-------
None
"""
# Open file for writing
f_dis = open(self.fn_path, "w")
# Item 0: heading
f_dis.write(f"{self.heading}\n")
# Item 1: NODES NLAY NJAG IVSD NPER ITMUNI LENUNI IDSYMRD
s = ""
for var in [
self.nodes,
self.nlay,
self.njag,
self.ivsd,
self.nper,
self.itmuni,
self.lenuni,
self.idsymrd,
]:
s += f"{var} "
f_dis.write(s + "\n")
# Item 2: LAYCBD
for k in range(self.nlay):
f_dis.write(f"{self.laycbd[k]:3d}")
f_dis.write("\n")
# Item 3: NODELAY
f_dis.write(self.nodelay.get_file_entry())
# Item 4: TOP
f_dis.write(self.top.get_file_entry())
# Item 5: BOT
f_dis.write(self.bot.get_file_entry())
# Item 6: AREA
f_dis.write(self.area.get_file_entry())
# Item 7: IAC
f_dis.write(self.iac.get_file_entry())
# Item 8: JA
f_dis.write(self.ja.get_file_entry())
# Item 9: IVC
if self.ivsd == 1:
f_dis.write(self.ivc.get_file_entry())
# Item 10a: CL1
if self.idsymrd == 1:
f_dis.write(self.cl1.get_file_entry())
# Item 10b: CL2
if self.idsymrd == 1:
f_dis.write(self.cl2.get_file_entry())
# Item 11: CL12
if self.idsymrd == 0:
f_dis.write(self.cl12.get_file_entry())
# Item 12: FAHL
f_dis.write(self.fahl.get_file_entry())
# Item 13: NPER, NSTP, TSMULT, Ss/tr
for t in range(self.nper):
f_dis.write(
f"{self.perlen[t]:14f}{self.nstp[t]:14d}{self.tsmult[t]:10f} "
)
if self.steady[t]:
f_dis.write(" SS\n")
else:
f_dis.write(" TR\n")
# Close and return
f_dis.close()
return
@staticmethod
def _ftype():
return "DISU"
@staticmethod
def _defaultunit():
return 11
def _get_neighboring_nodes(self):
"""
For each node, get node numbers for all neighbors.
Returns
-------
Jagged list of numpy arrays for each node.
Each array contains base-1 neighboring node indices.
"""
ja = self.ja.array
iac_sum = np.cumsum(self.iac.array)
ja_slices = np.asarray(
[
np.s_[iac_sum[i - 1] + 1 : x] if i > 0 else np.s_[1:x]
for i, x in enumerate(iac_sum)
]
) # note: this removes the diagonal - neighbors only
self._neighboring_nodes = [ja[sl] for sl in ja_slices]
return