import numpy as np
from .binarygrid_util import MfGrdFile
[docs]def get_structured_faceflows(
flowja, grb_file=None, ia=None, ja=None, verbose=False
):
"""
Get the face flows for the flow right face, flow front face, and
flow lower face from the MODFLOW 6 flowja flows. This method can
be useful for building face flow arrays for MT3DMS, MT3D-USGS, and
RT3D. This method only works for a structured MODFLOW 6 model.
Parameters
----------
flowja : ndarray
flowja array for a structured MODFLOW 6 model
grbfile : str
MODFLOW 6 binary grid file path
ia : list or ndarray
CRS row pointers. Only required if grb_file is not provided.
ja : list or ndarray
CRS column pointers. Only required if grb_file is not provided.
verbose: bool
Write information to standard output
Returns
-------
frf : ndarray
right face flows
fff : ndarray
front face flows
flf : ndarray
lower face flows
"""
if grb_file is not None:
grb = MfGrdFile(grb_file, verbose=verbose)
if grb.grid_type != "DIS":
raise ValueError(
"get_structured_faceflows method "
"is only for structured DIS grids"
)
ia, ja = grb.ia, grb.ja
else:
if ia is None or ja is None:
raise ValueError(
"ia and ja arrays must be specified if the MODFLOW 6"
"binary grid file name is not specified."
)
# flatten flowja, if necessary
if len(flowja.shape) > 0:
flowja = flowja.flatten()
# evaluate size of flowja relative to ja
__check_flowja_size(flowja, ja)
# create face flow arrays
shape = (grb.nlay, grb.nrow, grb.ncol)
frf = np.zeros(shape, dtype=float).flatten()
fff = np.zeros(shape, dtype=float).flatten()
flf = np.zeros(shape, dtype=float).flatten()
# fill flow terms
vmult = [-1.0, -1.0, -1.0]
flows = [frf, fff, flf]
for n in range(grb.nodes):
i0, i1 = ia[n] + 1, ia[n + 1]
ipos = 0
for j in range(i0, i1):
jcol = ja[j]
if jcol > n:
flows[ipos][n] = vmult[ipos] * flowja[j]
ipos += 1
# reshape flow terms
frf = frf.reshape(shape)
fff = fff.reshape(shape)
flf = flf.reshape(shape)
return frf, fff, flf
[docs]def get_residuals(
flowja, grb_file=None, ia=None, ja=None, shape=None, verbose=False
):
"""
Get the residual from the MODFLOW 6 flowja flows. The residual is stored
in the diagonal position of the flowja vector.
Parameters
----------
flowja : ndarray
flowja array for a structured MODFLOW 6 model
grbfile : str
MODFLOW 6 binary grid file path
ia : list or ndarray
CRS row pointers. Only required if grb_file is not provided.
ja : list or ndarray
CRS column pointers. Only required if grb_file is not provided.
shape : tuple
shape of returned residual. A flat array is returned if shape is None
and grbfile is None.
verbose: bool
Write information to standard output
Returns
-------
residual : ndarray
Residual for each cell
"""
if grb_file is not None:
grb = MfGrdFile(grb_file, verbose=verbose)
shape = grb.shape
ia, ja = grb.ia, grb.ja
else:
if ia is None or ja is None:
raise ValueError(
"ia and ja arrays must be specified if the MODFLOW 6 "
"binary grid file name is not specified."
)
# flatten flowja, if necessary
if len(flowja.shape) > 0:
flowja = flowja.flatten()
# evaluate size of flowja relative to ja
__check_flowja_size(flowja, ja)
# create residual
nodes = grb.nodes
residual = np.zeros(nodes, dtype=float)
# fill flow terms
for n in range(nodes):
i0, i1 = ia[n], ia[n + 1]
if i0 < i1:
residual[n] = flowja[i0]
else:
residual[n] = np.nan
# reshape residual terms
if shape is not None:
residual = residual.reshape(shape)
return residual
# internal
def __check_flowja_size(flowja, ja):
"""
Check the shape of flowja relative to ja.
"""
if flowja.shape != ja.shape:
raise ValueError(
f"size of flowja ({flowja.shape}) not equal to {ja.shape}"
)