Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix(check): Updated flopy's check to work with cellid -1 values #1885

Merged
merged 7 commits into from
Aug 6, 2023
Merged

fix(check): Updated flopy's check to work with cellid -1 values #1885

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c7a9dfe
fix(check): Updated check to work with cellid -1 values
scottrp Jul 27, 2023
dd5b7ac
fix(test): should be a complete grid
scottrp Jul 27, 2023
bced4eb
fix(model grid): Always resync mf2005 grids
scottrp Jul 27, 2023
80067bf
Revert "fix(model grid): Always resync mf2005 grids"
scottrp Jul 27, 2023
dcc62e6
fix(grid idomain): grid idomain now defaults to all 1's for MF6 only
scottrp Jul 28, 2023
12cfa82
fix(notebook): notebook contained out of bounds cellid
scottrp Jul 28, 2023
a9abc35
fix(modelgrid): always force modelgrid to resync until an idomain arr…
scottrp Jul 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions .docs/Notebooks/mf6_data_tutorial06.py
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,9 +102,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
View file
Open in desktop
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
View file
Open in desktop
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:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure I understand why there is a dependence of idomain on botm. Do we want to retain idomain=None if botm is not defined. Wouldn't it be better from the modelgrid creation to fail if botm has not been defined when creating a discretization instance?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Our current approach with modelgrid creation is to allow for the creation of a modelgrid without top and botm. This has always been allowed so that a user can make use of much of the modelgrid functionality without using top and botm. Also, a modelgrid is requested by flopy whenever grid data is loaded and the "check" option requested, including when top and botm are loaded. Failing to create modelgrid when botm is still in the process of loading would break our current approach. I am open to reconsidering this approach, but I thought that was beyond the scope of this PR, which is to fix a very specific problem.

Regarding the dependence of initializing the default idomain (of all 1's) on botm, to initialize the default idomain one needs to first know the shape of the model grid. The modelgrid class currently has only one variable, besides idomain, that defines the full shape of the grid, botm. The simplest, and hopefully least prone to error, implementation was therefore to get the shape of idomain from the botm array.

However, on further considering there is a potential problem with this. The model grid could get cached with the default idomain values when flopy loads botm and then not get updated when idomain is finally loaded. Going to add some code to fix this potential bug.

return idomain
else:
return np.ones_like(botm)
return idomain