import numpy as np
from ..modflow import Modflow
from .util_array import Util2d, Util3d
[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
----------
parent : flopy.modflow.Modflow
parent model
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 = idxl.min()
self.nplend = idxl.max()
self.npcbeg = idxc.min()
self.npcend = idxc.max()
self.nprbeg = idxr.min()
self.nprend = idxr.max()
# child grid dimensions
self.nlay = 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 + self.delrp[0 : self.npcbeg].sum()
self.yll = yllp + self.delcp[self.nprend + 1 :].sum()
return
[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] = self.delrp[j - 1] / 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] = self.delcp[i - 1] / 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 = 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 = pbotm[kp, ip, jp]
bot = 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 = 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 = topp[ip, jp]
btp = botp[kp, ip, jp]
if kp > 0:
tpp = botp[kp - 1, ip, jp]
tpc = topc[ic, jc]
btc = botc[kc, ic, jc]
if kc > 0:
tpc = botc[kc - 1, ic, jc]
if ihc == 0:
cl1 = 0.5 * (tpp - btp)
cl2 = 0.5 * (tpc - btc)
hwva = delrc[jc] * delcc[ic]
else:
if abs(idir) == 1:
cl1 = 0.5 * delrp[jp]
cl2 = 0.5 * delrc[jc]
hwva = delcc[ic]
elif abs(idir) == 2:
cl1 = 0.5 * delcp[ip]
cl2 = 0.5 * delcc[ic]
hwva = delrc[jc]
# connection distance
cd = None
if cdist:
if abs(idir) == 3:
cd = cl1 + cl2
else:
x1 = xc[ic, jc]
y1 = yc[ic, jc]
x2 = xp[ip, jp]
y2 = 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(angle)
if cdist:
exg.append(cd)
exglist.append(exg)
return exglist