import numpy as np
from ..discretization import StructuredGrid
from ..modflow import Modflow
from .util_array import Util2d, Util3d
[docs]class SimpleRegularGrid:
"""
Simple object for representing regular MODFLOW grid information.
Parameters
----------
nlay : int
number of layers
nrow : int
number of rows
ncol : int
number of columns
delr : ndarray
delr array
delc : ndarray
delc array
top : ndarray
top array (nrow, ncol)
botm : ndarray
botm array (nlay, nrow, ncol)
idomain : ndarray
idomain array (nlay, nrow, ncol)
xorigin : float
x location of grid lower left corner
yorigin : float
y location of grid lower left corner
"""
def __init__(
self,
nlay,
nrow,
ncol,
delr,
delc,
top,
botm,
idomain,
xorigin,
yorigin,
):
# enforce compliance
assert delr.shape == (ncol,)
assert delc.shape == (nrow,)
assert top.shape == (nrow, ncol)
assert botm.shape == (nlay, nrow, ncol)
assert idomain.shape == (nlay, nrow, ncol)
self.nlay = nlay
self.nrow = nrow
self.ncol = ncol
self.delr = delr
self.delc = delc
self.top = top
self.botm = botm
self.idomain = idomain
self.xorigin = xorigin
self.yorigin = yorigin
return
@property
def modelgrid(self):
mg = StructuredGrid(
delc=self.delc,
delr=self.delr,
top=self.top,
botm=self.botm,
idomain=self.idomain,
xoff=self.xorigin,
yoff=self.yorigin,
)
return mg
[docs] def get_gridprops_dis6(self):
gridprops = {
"xorigin": self.xorigin,
"yorigin": self.yorigin,
"nlay": self.nlay,
"nrow": self.nrow,
"ncol": self.ncol,
"delr": self.delr,
"delc": self.delc,
"top": self.top,
"botm": self.botm,
"idomain": self.idomain,
}
return gridprops
[docs]class Lgr:
def __init__(
self,
nlayp,
nrowp,
ncolp,
delrp,
delcp,
topp,
botmp,
idomainp,
ncpp=3,
ncppl=1,
xllp=0.0,
yllp=0.0,
):
"""
Parameters
----------
nlayp : int
parent layers
nrowp : int
parent number of rows
ncolp : int
parent number of columns
delrp : ndarray
parent delr array
delcp : ndarray
parent delc array
topp : ndarray
parent top array (nrowp, ncolp)
botmp : ndarray
parent botm array (nlayp, nrowp, ncolp)
idomainp : ndarray
parent idomain array used to create the child grid. Ones indicate
a parent cell and zeros indicate a child cell. The domain of the
child grid will span a rectangular region that spans all idomain
cells with a value of zero. idomain must be of shape
(nlayp, nrowp, ncolp)
ncpp : int
number of child cells along the face of a parent cell
ncppl : list of ints
number of child layers per parent layer
xllp : float
x location of parent grid lower left corner
yllp : float
y location of parent grid lower left corner
"""
# parent grid properties
self.nlayp = nlayp
self.nrowp = nrowp
self.ncolp = ncolp
m = Modflow()
self.delrp = Util2d(m, (ncolp,), np.float32, delrp, "delrp").array
self.delcp = Util2d(m, (nrowp,), np.float32, delcp, "delcp").array
self.topp = Util2d(m, (nrowp, ncolp), np.float32, topp, "topp").array
self.botmp = Util3d(
m, (nlayp, nrowp, ncolp), np.float32, botmp, "botmp"
).array
# idomain
assert idomainp.shape == (nlayp, nrowp, ncolp)
self.idomain = idomainp
idxl, idxr, idxc = np.where(idomainp == 0)
assert idxl.shape[0] > 1, "no zero values found in idomain"
# # child cells per parent and child cells per parent layer
self.ncpp = ncpp
self.ncppl = Util2d(m, (nlayp,), np.int32, ncppl, "ncppl").array
# calculate ibcl which is the bottom child layer (one based) in each
# parent layer
self.ibcl = np.zeros(self.nlayp, dtype=int)
self.ibcl[0] = self.ncppl[0]
for k in range(1, self.nlayp):
if self.ncppl[k] > 0:
self.ibcl[k] = self.ibcl[k - 1] + self.ncppl[k]
# parent lower left
self.xllp = xllp
self.yllp = yllp
# child grid properties
self.nplbeg = int(idxl.min())
self.nplend = int(idxl.max())
self.npcbeg = int(idxc.min())
self.npcend = int(idxc.max())
self.nprbeg = int(idxr.min())
self.nprend = int(idxr.max())
# child grid dimensions
self.nlay = int(self.ncppl.sum())
self.nrow = (self.nprend - self.nprbeg + 1) * ncpp
self.ncol = (self.npcend - self.npcbeg + 1) * ncpp
# assign child properties
self.delr, self.delc = self.get_delr_delc()
self.top, self.botm = self.get_top_botm()
self.xll = xllp + float(self.delrp[0 : self.npcbeg].sum())
self.yll = yllp + float(self.delcp[self.nprend + 1 :].sum())
[docs] def get_shape(self):
"""
Return the shape of the child grid
Returns
-------
(nlay, nrow, ncol) : tuple
shape of the child grid
"""
return self.nlay, self.nrow, self.ncol
[docs] def get_lower_left(self):
"""
Return the lower left corner of the child grid
Returns
-------
(xll, yll) : tuple
location of lower left corner of the child grid
"""
return self.xll, self.yll
[docs] def get_delr_delc(self):
# create the delr and delc arrays for this child grid
delr = np.zeros((self.ncol), dtype=float)
delc = np.zeros((self.nrow), dtype=float)
jstart = 0
jend = self.ncpp
for j in range(self.npcbeg, self.npcend + 1):
delr[jstart:jend] = float(self.delrp[j]) / self.ncpp
jstart = jend
jend = jstart + self.ncpp
istart = 0
iend = self.ncpp
for i in range(self.nprbeg, self.nprend + 1):
delc[istart:iend] = float(self.delcp[i]) / self.ncpp
istart = iend
iend = istart + self.ncpp
return delr, delc
[docs] def get_top_botm(self):
bt = self.botmp
tp = self.topp
shp = tp.shape
tp = tp.reshape(1, shp[0], shp[1])
pbotm = np.vstack((tp, bt))
botm = np.zeros((self.nlay + 1, self.nrow, self.ncol), dtype=float)
for ip in range(self.nprbeg, self.nprend + 1):
for jp in range(self.npcbeg, self.npcend + 1):
top = float(pbotm[0, ip, jp])
icrowstart = (ip - self.nprbeg) * self.ncpp
icrowend = icrowstart + self.ncpp
iccolstart = (jp - self.npcbeg) * self.ncpp
iccolend = iccolstart + self.ncpp
botm[0, icrowstart:icrowend, iccolstart:iccolend] = top
kc = 1
for kp in range(self.nplbeg, self.nplend + 1):
top = float(pbotm[kp, ip, jp])
bot = float(pbotm[kp + 1, ip, jp])
dz = (top - bot) / self.ncppl[kp]
for _ in range(self.ncppl[kp]):
botm[kc, icrowstart:icrowend, iccolstart:iccolend] = (
botm[
kc - 1,
icrowstart:icrowend,
iccolstart:iccolend,
]
- dz
)
kc += 1
return botm[0], botm[1:]
[docs] def get_replicated_parent_array(self, parent_array):
"""
Get a two-dimensional array the size of the child grid that has values
replicated from the provided parent array.
Parameters
----------
parent_array : ndarray
A two-dimensional array that is the size of the parent model rows
and columns.
Returns
-------
child_array : ndarray
A two-dimensional array that is the size of the child model rows
and columns
"""
assert parent_array.shape == (self.nrowp, self.ncolp)
child_array = np.empty(
(self.nrow, self.ncol), dtype=parent_array.dtype
)
for ip in range(self.nprbeg, self.nprend + 1):
for jp in range(self.npcbeg, self.npcend + 1):
icrowstart = (ip - self.nprbeg) * self.ncpp
icrowend = icrowstart + self.ncpp
iccolstart = (jp - self.npcbeg) * self.ncpp
iccolend = iccolstart + self.ncpp
value = int(parent_array[ip, jp])
child_array[icrowstart:icrowend, iccolstart:iccolend] = value
return child_array
[docs] def get_idomain(self):
"""
Return the idomain array for the child model. This will normally
be all ones unless the idomain array for the parent model is
non-rectangular and irregularly shaped. Then, parts of the child
model will have idomain zero cells.
Returns
-------
idomain : ndarray
idomain array for the child model
"""
idomain = np.ones((self.nlay, self.nrow, self.ncol), dtype=int)
for kc in range(self.nlay):
for ic in range(self.nrow):
for jc in range(self.ncol):
kp, ip, jp = self.get_parent_indices(kc, ic, jc)
if self.idomain[kp, ip, jp] == 1:
idomain[kc, ic, jc] = 0
return idomain
[docs] def get_parent_indices(self, kc, ic, jc):
"""
Method returns the parent cell indices for this child.
The returned indices are in zero-based indexing.
"""
ip = self.nprbeg + int(ic / self.ncpp)
jp = self.npcbeg + int(jc / self.ncpp)
kp = 0
kcstart = 0
for k in range(self.nplbeg, self.nplend + 1):
kcend = kcstart + self.ncppl[k] - 1
if kcstart <= kc <= kcend:
kp = k
break
kcstart = kcend + 1
return kp, ip, jp
[docs] def get_parent_connections(self, kc, ic, jc):
"""
Return a list of parent cell indices that are connected to child
cell kc, ic, jc.
"""
assert 0 <= kc < self.nlay, "layer must be >= 0 and < child nlay"
assert 0 <= ic < self.nrow, "layer must be >= 0 and < child nrow"
assert 0 <= jc < self.ncol, "layer must be >= 0 and < child ncol"
parentlist = []
(kp, ip, jp) = self.get_parent_indices(kc, ic, jc)
# parent cell to left
if jc % self.ncpp == 0:
if jp - 1 >= 0:
if self.idomain[kp, ip, jp - 1] != 0:
parentlist.append(((kp, ip, jp - 1), -1))
# parent cell to right
if (jc + 1) % self.ncpp == 0:
if jp + 1 < self.ncolp:
if self.idomain[kp, ip, jp + 1] != 0:
parentlist.append(((kp, ip, jp + 1), 1))
# parent cell to back
if ic % self.ncpp == 0:
if ip - 1 >= 0:
if self.idomain[kp, ip - 1, jp] != 0:
parentlist.append(((kp, ip - 1, jp), 2))
# parent cell to front
if (ic + 1) % self.ncpp == 0:
if ip + 1 < self.nrowp:
if self.idomain[kp, ip + 1, jp] != 0:
parentlist.append(((kp, ip + 1, jp), -2))
# parent cell to top is not possible
# parent cell to bottom
if kc + 1 == self.ibcl[kp]:
if kp + 1 < self.nlayp:
if self.idomain[kp + 1, ip, jp] != 0:
parentlist.append(((kp + 1, ip, jp), -3))
return parentlist
[docs] def get_exchange_data(self, angldegx=False, cdist=False):
"""
Get the list of parent/child connections
<cellidm1> <cellidm2> <ihc> <cl1> <cl2> <hwva> <angledegx>
Returns
-------
exglist : list
list of connections between parent and child
"""
exglist = []
nlayc = self.nlay
nrowc = self.nrow
ncolc = self.ncol
delrc = self.delr
delcc = self.delc
delrp = self.delrp
delcp = self.delcp
topp = self.topp
botp = self.botmp
topc = self.top
botc = self.botm
if cdist:
# child xy meshgrid
xc = np.add.accumulate(delrc) - 0.5 * delrc
Ly = np.add.reduce(delcc)
yc = Ly - (np.add.accumulate(delcc) - 0.5 * delcc)
xc += self.xll
yc += self.yll
xc, yc = np.meshgrid(xc, yc)
# parent xy meshgrid
xp = np.add.accumulate(delrp) - 0.5 * delrp
Ly = np.add.reduce(delcp)
yp = Ly - (np.add.accumulate(delcp) - 0.5 * delcp)
xc += self.xllp
yc += self.yllp
xp, yp = np.meshgrid(xp, yp)
cidomain = self.get_idomain()
for kc in range(nlayc):
for ic in range(nrowc):
for jc in range(ncolc):
plist = self.get_parent_connections(kc, ic, jc)
for (kp, ip, jp), idir in plist:
if cidomain[kc, ic, jc] == 0:
continue
# horizontal or vertical connection
# 1 if a child cell horizontally connected to a parent
# cell
# 2 if more than one child cells horizontally connected
# to parent cell
# 0 if a vertical connection
ihc = 1
if self.ncppl[kp] > 1:
ihc = 2
if abs(idir) == 3:
ihc = 0
# angldegx
angle = None
if angldegx:
angle = 180.0 # -x, west
if idir == 2:
angle = 270.0 # -y, south
elif idir == -1:
angle = 0.0 # +x, east
elif idir == -2:
angle = 90.0 # +y, north
# vertical connection
cl1 = None
cl2 = None
hwva = None
tpp = float(topp[ip, jp])
btp = float(botp[kp, ip, jp])
if kp > 0:
tpp = float(botp[kp - 1, ip, jp])
tpc = float(topc[ic, jc])
btc = float(botc[kc, ic, jc])
if kc > 0:
tpc = float(botc[kc - 1, ic, jc])
if ihc == 0:
cl1 = 0.5 * (tpp - btp)
cl2 = 0.5 * (tpc - btc)
hwva = float(delrc[jc]) * float(delcc[ic])
else:
if abs(idir) == 1:
cl1 = 0.5 * float(delrp[jp])
cl2 = 0.5 * float(delrc[jc])
hwva = float(delcc[ic])
elif abs(idir) == 2:
cl1 = 0.5 * float(delcp[ip])
cl2 = 0.5 * float(delcc[ic])
hwva = delrc[jc]
# connection distance
cd = None
if cdist:
if abs(idir) == 3:
cd = cl1 + cl2
else:
x1 = float(xc[ic, jc])
y1 = float(yc[ic, jc])
x2 = float(xp[ip, jp])
y2 = float(yp[ip, jp])
cd = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
exg = [
(kp, ip, jp),
(kc, ic, jc),
ihc,
cl1,
cl2,
hwva,
]
if angldegx:
exg.append(float(angle))
if cdist:
exg.append(float(cd))
exglist.append(exg)
return exglist
@property
def parent(self):
"""
Return a SimpleRegularGrid object for the parent model
Returns
-------
simple_regular_grid : SimpleRegularGrid
simple grid object containing grid information for the parent
"""
simple_regular_grid = SimpleRegularGrid(
self.nlayp,
self.nrowp,
self.ncolp,
self.delrp,
self.delcp,
self.topp,
self.botmp,
self.idomain,
self.xllp,
self.yllp,
)
return simple_regular_grid
@property
def child(self):
"""
Return a SimpleRegularGrid object for the child model
Returns
-------
simple_regular_grid : SimpleRegularGrid
simple grid object containing grid information for the child
"""
delrc, delcc = self.get_delr_delc()
idomainc = self.get_idomain() # child idomain
topc = self.top
botmc = self.botm
child_dis_shp = self.get_shape()
nlayc = child_dis_shp[0]
nrowc = child_dis_shp[1]
ncolc = child_dis_shp[2]
xorigin = self.xll
yorigin = self.yll
simple_regular_grid = SimpleRegularGrid(
nlayc,
nrowc,
ncolc,
delrc,
delcc,
topc,
botmc,
idomainc,
xorigin,
yorigin,
)
return simple_regular_grid