Skip to content

Commit

Permalink
fix(check): Updated flopy's check to work with cellid -1 values (#1885)
Browse files Browse the repository at this point in the history
* fix(test): should be a complete grid
* fix(model grid): Always resync mf2005 grids
* Revert "fix(model grid): Always resync mf2005 grids"
* fix(grid idomain): grid idomain now defaults to all 1's for MF6 only
* fix(notebook): notebook contained out of bounds cellid
* fix(modelgrid): always force modelgrid to resync until an idomain array is provided

---------

Co-authored-by: scottrp <[email protected]>
  • Loading branch information
spaulins-usgs and scottrp authored Aug 6, 2023
1 parent 8084b67 commit 68f0c3b
Show file tree
Hide file tree
Showing 5 changed files with 71 additions and 17 deletions.
9 changes: 5 additions & 4 deletions .docs/Notebooks/mf6_data_tutorial06.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
# ---
# jupyter:
# jupytext:
# notebook_metadata_filter: metadata
# text_representation:
# extension: .py
# format_name: light
# format_version: "1.5"
# jupytext_version: 1.5.1
# format_version: '1.5'
# jupytext_version: 1.14.4
# kernelspec:
# display_name: Python 3
# display_name: Python 3 (ipykernel)
# language: python
# name: python3
# metadata:
Expand Down Expand Up @@ -159,7 +160,7 @@
# Note that the cellid information (layer, row, column) is encapsulated in
# a tuple.

stress_period_data = [((1, 10, 10), 100.0), ((1, 10, 11), 105.0)]
stress_period_data = [((1, 8, 8), 100.0), ((1, 9, 9), 105.0)]
# build chd package
chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(
gwf,
Expand Down
4 changes: 3 additions & 1 deletion autotest/regression/test_mf6.py
Original file line number Diff line number Diff line change
Expand Up @@ -677,6 +677,7 @@ def test_np001(function_tmpdir, example_data_path):
sim.delete_output_files()

# test error checking
sim.simulation_data.verify_data = False
drn_package = ModflowGwfdrn(
model,
print_input=True,
Expand All @@ -697,6 +698,7 @@ def test_np001(function_tmpdir, example_data_path):
k=100001.0,
k33=1e-12,
)
sim.simulation_data.verify_data = True
chk = sim.check()
summary = ".".join(chk[0].summary_array.desc)
assert "drn_1 package: invalid BC index" in summary
Expand Down Expand Up @@ -3083,7 +3085,7 @@ def test028_create_tests_sfr(function_tmpdir, example_data_path):
delc=5000.0,
top=top,
botm=botm,
idomain=idomain,
#idomain=idomain,
filename=f"{model_name}.dis",
)
strt = testutils.read_std_array(os.path.join(pth, "strt.txt"), "float")
Expand Down
3 changes: 0 additions & 3 deletions flopy/discretization/vertexgrid.py
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,6 @@ def __init__(
self._vertices = vertices
self._cell1d = cell1d
self._cell2d = cell2d
self._top = top
self._botm = botm
self._idomain = idomain
if botm is None:
self._nlay = nlay
self._ncpl = ncpl
Expand Down
32 changes: 31 additions & 1 deletion flopy/mf6/data/mfdatalist.py
Original file line number Diff line number Diff line change
Expand Up @@ -505,6 +505,36 @@ def _check_valid_cellids(self):
idomain_val = idomain
# cellid should be within the model grid
for idx, cellid_part in enumerate(record[index]):
if cellid_part == -1:
# cellid not defined, all values should
# be -1
match = all(
elem == record[index][0]
for elem in record[index]
)
if not match:
message = (
f"Invalid cellid {record[index]}"
)
(
type_,
value_,
traceback_,
) = sys.exc_info()
raise MFDataException(
self.structure.get_model(),
self.structure.get_package(),
self.structure.path,
"storing data",
self.structure.name,
inspect.stack()[0][3],
type_,
value_,
traceback_,
message,
self._simulation_data.debug,
)
continue
if (
model_shape[idx] <= cellid_part
or cellid_part < 0
Expand All @@ -530,7 +560,7 @@ def _check_valid_cellids(self):
)
idomain_val = idomain_val[cellid_part]
# cellid should be at an active cell
if idomain_val < 1:
if record[index][0] != -1 and idomain_val < 1:
message = (
"Cellid {} is outside of the "
"active model grid"
Expand Down
40 changes: 32 additions & 8 deletions flopy/mf6/mfmodel.py
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,7 @@ def modelgrid(self):
model.
"""

force_resync = False
if not self._mg_resync:
return self._modelgrid
if self.get_grid_type() == DiscretizationType.DIS:
Expand All @@ -370,12 +370,17 @@ def modelgrid(self):
angrot=self._modelgrid.angrot,
)
else:
botm = dis.botm.array
idomain = dis.idomain.array
if idomain is None:
force_resync = True
idomain = self._resolve_idomain(idomain, botm)
self._modelgrid = StructuredGrid(
delc=dis.delc.array,
delr=dis.delr.array,
top=dis.top.array,
botm=dis.botm.array,
idomain=dis.idomain.array,
botm=botm,
idomain=idomain,
lenuni=dis.length_units.array,
crs=self._modelgrid.crs,
xoff=self._modelgrid.xoffset,
Expand Down Expand Up @@ -403,12 +408,17 @@ def modelgrid(self):
angrot=self._modelgrid.angrot,
)
else:
botm = dis.botm.array
idomain = dis.idomain.array
if idomain is None:
force_resync = True
idomain = self._resolve_idomain(idomain, botm)
self._modelgrid = VertexGrid(
vertices=dis.vertices.array,
cell2d=dis.cell2d.array,
top=dis.top.array,
botm=dis.botm.array,
idomain=dis.idomain.array,
botm=botm,
idomain=idomain,
lenuni=dis.length_units.array,
crs=self._modelgrid.crs,
xoff=self._modelgrid.xoffset,
Expand Down Expand Up @@ -498,12 +508,17 @@ def modelgrid(self):
angrot=self._modelgrid.angrot,
)
else:
botm = dis.botm.array
idomain = dis.idomain.array
if idomain is None:
force_resync = True
idomain = self._resolve_idomain(idomain, botm)
self._modelgrid = VertexGrid(
vertices=dis.vertices.array,
cell1d=dis.cell1d.array,
top=dis.top.array,
botm=dis.botm.array,
idomain=dis.idomain.array,
botm=botm,
idomain=idomain,
lenuni=dis.length_units.array,
crs=self._modelgrid.crs,
xoff=self._modelgrid.xoffset,
Expand Down Expand Up @@ -546,7 +561,7 @@ def modelgrid(self):
angrot,
self._modelgrid.crs,
)
self._mg_resync = not self._modelgrid.is_complete
self._mg_resync = not self._modelgrid.is_complete or force_resync
return self._modelgrid

@property
Expand Down Expand Up @@ -1937,3 +1952,12 @@ def plot(self, SelPackList=None, **kwargs):
)

return axes

@staticmethod
def _resolve_idomain(idomain, botm):
if idomain is None:
if botm is None:
return idomain
else:
return np.ones_like(botm)
return idomain

0 comments on commit 68f0c3b

Please sign in to comment.