modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
507 stars 307 forks source link

bug: `modelgrid` coord info dropped when using a DISV grid #2283

Closed cneyens closed 1 month ago

cneyens commented 1 month ago

Describe the bug When obtaining the modelgrid object from a model with a DISV package, the coord info present in the DISV object is not passed properly to the modelgrid object. For models with DIS, this works as expected; I haven't checked for DISU grids.

To Reproduce Example model with a triangular grid:

Example code ```py import os import numpy as np import flopy from flopy.utils.triangle import Triangle # create triangular grid triangle_ws = 'triangle' if not os.path.exists(triangle_ws): os.mkdir(triangle_ws) active_area = [(0, 0), (0, 1000), (1000, 1000), (1000, 0)] tri = Triangle(model_ws=triangle_ws, angle=30) tri.add_polygon(active_area) tri.add_region((1, 1), maximum_area=50**2) tri.build() # build vertex grid object vgrid = flopy.discretization.VertexGrid( vertices=tri.get_vertices(), cell2d=tri.get_cell2d(), xoff=199000, yoff=215500, crs=31370, angrot=30, ) # coord info is set (also correct when using vgrid.set_coord_info() print(vgrid) # create MODFLOW 6 model sim = flopy.mf6.MFSimulation(sim_name='prj-test', sim_ws='model') tdis = flopy.mf6.ModflowTdis(sim) ims = flopy.mf6.ModflowIms(sim) gwf = flopy.mf6.ModflowGwf(sim, modelname='gwf') disv = flopy.mf6.ModflowGwfdisv( gwf, xorigin=vgrid.xoffset, yorigin=vgrid.yoffset, angrot=vgrid.angrot, # no CRS info can be set in DISV nlay=1, top=0.0, botm=-10.0, ncpl=vgrid.ncpl, nvert=vgrid.nvert, cell2d=vgrid.cell2d, vertices=tri.get_vertices() # this is not stored in the Vertex grid object? ) # coord info is dropped # also the case when loading model from disk print(gwf.modelgrid) # disv object has coord info set properly though print([disv.xorigin.get_data(), disv.yorigin.get_data(), disv.angrot.get_data()]) ```

Expected behavior I expect the modelgrid object to get the coord info from the DISV object when calling it through gwf.modelgrid, as is done for a DIS model. Pretty easy to get around this using set_coord_info() based on the DISV properties, but perhaps there's more to it.

Desktop (please complete the following information):

martclanor commented 1 month ago

This seems to be caused by this block of code in flopy/mf6/mfmodel.py:575 (https://github.com/modflowpy/flopy/blob/develop/flopy/mf6/mfmodel.py#L575)

if self.get_grid_type() != DiscretizationType.DISV:
    # get coordinate data from dis file
    xorig = dis.xorigin.get_data()
    yorig = dis.yorigin.get_data()
    angrot = dis.angrot.get_data()
else:
    xorig = self._modelgrid.xoffset
    yorig = self._modelgrid.yoffset
    angrot = self._modelgrid.angrot

I am tempted to just remove this exception to the DISV case, but I'm not sure yet what it is intended for.

langevin-usgs commented 1 month ago

I can't think of a reason why this was handled this way, except that xorigin, yorigin, and angrot weren't originally included with DISV and DISU.

image