feat(grid): add thickness method for all grid types

[jdhughes-dev]

Includes option for passing array to calculate saturated thickness.

[jdhughes-dev]

@jlarsen-usgs Added the proposed grid.thick() method for Mike Fienen.

It also seems the grid.top_botm property was only appropriate for structured grids.

Let me know what you think of the proposed changes. If it looks good I will edit the pull request and deprecate the thickness property on the ModflowDis class (we should probably do the same for the zcentroids property) and point users to the grid class.

[codecov]

Codecov Report

Merging #1138 (23edc80) into develop (fbba027) will decrease coverage by 0.123%.
The diff coverage is 82.242%.

❗ Current head 23edc80 differs from pull request most recent head d341ede. Consider uploading reports for the commit d341ede to get more accurate results

@@              Coverage Diff              @@
##           develop     #1138       +/-   ##
=============================================
- Coverage   71.729%   71.605%   -0.124%     
=============================================
  Files          225       225               
  Lines        51919     51967       +48     
=============================================
- Hits         37241     37211       -30     
- Misses       14678     14756       +78     

Impacted Files	Coverage Δ
flopy/mf6/modflow/mfgwf.py	55.555% <ø> (ø)
flopy/modflow/mfdisu.py	80.981% <0.000%> (ø)
flopy/plot/crosssection.py	60.225% <ø> (+4.513%)	⬆️
flopy/utils/postprocessing.py	75.510% <42.857%> (-7.726%)	⬇️
flopy/modflow/mfdis.py	86.783% <50.000%> (+0.033%)	⬆️
flopy/discretization/unstructuredgrid.py	70.491% <57.142%> (+5.158%)	⬆️
flopy/discretization/vertexgrid.py	77.083% <76.470%> (+3.989%)	⬆️
flopy/discretization/structuredgrid.py	87.466% <88.235%> (+3.428%)	⬆️
flopy/discretization/grid.py	74.378% <90.000%> (+0.388%)	⬆️
flopy/plot/map.py	70.748% <94.736%> (ø)
... and 19 more

[jlarsen-usgs]

Hi Joe,

I've added some general comments and changes that will help remove the if/else statements from the code. The main being implementing a top_botm property method in the UnstructuredGrid class that overrides the base Grid class.

[jdhughes-dev]

@jlarsen-usgs I implemented your suggestions. I modified them a bit because the unstructured was returning an array for the thickness instead of a vector. So thick return arrays with (nlay, nrow. ncol), (nlay, ncpl), and (nodes) shapes for structured, vertex, and unstructured grids.

Also made thick a property and added a separate saturated_thick method for calculating the saturated thickness.

[jdhughes-dev]

@aleaf I was working on a generic thick property for the grid class based on a conversation I had with @mnfienen yesterday. I also added a generic saturated_thick() method.

When I was looking into replacing get_saturated_thickness() in t042 with the new grid class saturated_thick() method. I noticed that the current get_saturated_thickness() method calculates a value of 1.6 for cell (2,1,1) which has top and bottom elevations of 2 and 1, respectively, and a head of 2.6. This isn't what I would call the saturated thickness which to me is 1.0. Instead I would say get_saturated_thickness() is calculating a pressure head for this case.

I found what looks to be a bug in the logic in get_daturated_thickness(). I have modified

unconf_thickness = np.where((hds - botm) > top, top, hds - botm)

to

unconf_thickness = np.where(hds > top, top - botm, hds - botm)

I couldn't understand why a saturated thickness would be compared to the top elevation.

Also added a check for heads less than the bottom of the cell (for MODFLOW-NWT). Added the following above the previous statement.

hds = np.where(hds < botm, botm, hds)  # for NWT when hds < botm

[jdhughes-dev]

@jlarsen-usgs You should also look at the t050 and t072 autotests (I think you are the one who worked on the parts of the autotests that I modified). I had to change some of the values in the assertions. The specific discharge values are different because I replace the postprocessing,get_saturated_thickness() method with the grid.saturated_thick() method. For cases where the head is above the top of the cell the postprocessing,get_saturated_thickness() method is calculating the pressure head not the aquifer thickness. So the fact that the vectors are changes does not surprise me but you should look at it to confirm this.

[jlarsen-usgs]

@jdhughes-usgs

I looked through the new changes for calculating saturated thickness and the associated tests and they make sense to me. The changes in vector values are what I would expect from changing from potentiometric surface to layer thickness when heads are higher than a given model cell's top elevation.