spatialreference not loading correctly

[jtwhite79]

I know, I know, modelgrid and whatnot, but a lot of us still rely on SpatialReference. This code shows the problem:

import flopy
m = flopy.modflow.Modflow("test",xll=12345,yll=12345)
flopy.modflow.ModflowDis(m,10,10,10)
m.sr.xll = 12345
m.sr.yll = 12345
print(m.sr)
print(m.modelgrid)
m.write_input()
mm = flopy.modflow.Modflow.load("test.nam")
print(mm.sr)
print(mm.modelgrid)

The loaded model (mm) does not have the right sr information.

So Im sure someone will tell me to get with the times, which I have tried. I tried to use modelgrid, but unlike sr, which caches the xcentergrid and ycentergrid arrays after first call, modelgrid.xcellcenters is dynamically refilling the arrays with every call! I guess I get the arrays I need once and then process with them, but I like the encapsulation of caching them for users - this breaks the idea of "just find and replace sr with modelgrid". Anyway, can someone please fix sr?

[christianlangevin]

@jtwhite79, sounds like @jlarsen-usgs will take a look. Those x and y cell centers should be cached. That was something that @spaulins-usgs added a while back, I think.

[jtwhite79]

ok - I could be mistaken about the caching - I just noticed my computer was about to melt when I switched to modelgrid - the following ctrl+c showed a traceback in the vertices calculations...

[spaulins-usgs]

@jtwhite79, the grid class caches the results from your last xcellcenters call as long as you don't make changes to the grid. Calls like set_coord_info clears the cache. Also, by default xcellcenters returns a copy of the cached data instead of a reference to that cached data. It was designed this way to prevent users from accidentally overwriting the cache. Grid has a hidden _copy_cache attribute that can be turned off in order to get a reference to the data.

Also when you call xcellcenters multiple times make sure each call is working with the same grid object. The model object does not cache a model grid, but constructs a new one for you every time you call .modelgrid. This is to guarantee that the a newly requested grid object is sync with the discritization package. Since the model does not cache the grid this code:

for n in range(0, 100):
    m.modelgrid.xcellcenters

does not cache the xcellcenters results and may take significantly longer to run than the following code which does cache xcellcenters results:

mg = m.modelgrid
for n in range(0, 100):
    mg.xcellcenters

[jtwhite79]

thanks @spaulins-usgs - that first option is more or less what I was trying to do.

I can see a few use cases where the discret might change relative to the modelgrid, but (at least in my narrow use of flopy) seems like the more common case would be that the discret and modelgrid are static from instantiation, right? Im afraid the dynamic modelgrid creation is going to crush performance unless users remember this trick (at least that is what happened to me...).

[spaulins-usgs]

I agree it is much more common that discret and modelgrid are static from instantiation. Unfortunately if I assume that they are static and always provide the same cached model grid it will either produce incorrect results or crash in the few cases where they are not static.

One solution is to have the model class check discret against a cached model grid and provide either the cached modelgrid or a new one based on if anything has changed. Though doing this may just replace one performance issue with another since verifying the consistency of potentially large arrays like botm and idomain might be slow.

Another solution might be to provide a warning message when a user requests a modelgrid more than once from the same method. The warning message could be something like:

"WARNING: Property .modelgrid called multiple times from the same method. Each time .modelgrid is called it creates a new grid object. If this is not what you intended to do and you are seeing performance issues try working with the model grid as a local variable 'mg = model.modegrid'".

I am leaning toward the second option though I am open to suggestions.

[christianlangevin]

@spaulins-usgs, isn't there any way to monitor for changes that would require modelgrid to be recreated? And I don't mean checking the array values against each other, because I agree that would be slow, I mean having an attribute set for dis (or packages in general) that are reset whenever something changes. Then, if the DIS package changes or if something changes in modelgrid, then it knows to rebuild it. Seems like we had something like that in spatial reference.

[christianlangevin]

BTW, @jtwhite79, I think this issue should be fixed so I'm closing. But would be good to have on your input on all this so we can keep moving forward...