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

Add saturated_thick method to grid class to calculate saturated
thickness. Can pass a mask value(s) to filter out cells. Added
deprecation warning to ModflowDis for thickness property and
postprocessing.get_saturated_thickness(). Revised notebooks and
autotests to use grid.thick property and grid.saturated_thick method.
Fix bug in postprocessing.get_saturated_thickness() method.

autotest/t042_test.py

@@ -116,10 +116,10 @@ def test_get_sat_thickness_gradients():
     dz = np.array([np.nan, -0.9])
     assert np.nansum(np.abs(dh / dz - grad[:, 1, 0])) < 1e-6
 
-    sat_thick = get_saturated_thickness(hds, m, nodata)
+    sat_thick = m.modelgrid.saturated_thick(hds, mask=nodata)
     assert (
-        np.abs(np.sum(sat_thick[:, 1, 1] - np.array([0.2, 1.0, 1.6]))) < 1e-6
-    )
+        np.abs(np.sum(sat_thick[:, 1, 1] - np.array([0.2, 1.0, 1.0]))) < 1e-6
+    ), "failed saturated thickness comparison (grid.thick())"
 
 
 if __name__ == "__main__":

autotest/t050_test.py

@@ -468,22 +468,26 @@ def test_vtk_vector():
     # assert(totalbytes1==942413)
     # nlines1 = count_lines_in_file(filetocheck, binary=True)
     # assert(nlines1==3824)
-    assert os.path.exists(filetocheck)
+    assert os.path.exists(filetocheck), "file (0) does not exist: {}".format(
+        filetocheck
+    )
 
     # with values directly given at vertices
     q = pp.get_specific_discharge(vectors, m, head, position="vertices")
     nancount = np.count_nonzero(np.isnan(q[0]))
-    assert nancount == 472
+    assert nancount == 472, "nancount != 472 ({})".format(nancount)
     overall = np.nansum(q[0]) + np.nansum(q[1]) + np.nansum(q[2])
-    assert np.allclose(overall, -45.38671967357735)
+    assert np.allclose(
+        overall, -15.849639024891047
+    ), "vertices overall = {}".format(overall)
     output_dir = os.path.join(cpth, "freyberg_vector")
     filenametocheck = "discharge_verts.vtu"
     vtk.export_vector(m, q, output_dir, "discharge_verts")
     filetocheck = os.path.join(output_dir, filenametocheck)
     # totalbytes2 = os.path.getsize(filetocheck)
     # assert(totalbytes2==1990047)
     nlines2 = count_lines_in_file(filetocheck)
-    assert nlines2 == 10598
+    assert nlines2 == 10598, "nlines != 10598 ({})".format(nlines2)
 
     # with values directly given at vertices and binary
     vtk.export_vector(
@@ -494,7 +498,9 @@ def test_vtk_vector():
     # assert(totalbytes3==891486)
     # nlines3 = count_lines_in_file(filetocheck, binary=True)
     # assert(nlines3==3012)
-    assert os.path.exists(filetocheck)
+    assert os.path.exists(filetocheck), "file (1) does not exist: {}".format(
+        filetocheck
+    )
 
     return
 

autotest/t072_test_spedis.py

@@ -450,11 +450,11 @@ def test_specific_discharge_comprehensive():
     for col in ax.collections:
         if not isinstance(col, Quiver):
             raise AssertionError("Unexpected collection type")
-    assert np.sum(quiver.Umask) == 2
+    assert np.sum(quiver.Umask) == 4
     pos = np.sum(quiver.X) + np.sum(quiver.Y)
     assert np.allclose(pos, 1600.0)
     val = np.sum(quiver.U) + np.sum(quiver.V)
-    assert np.allclose(val, 10.11359908150753)
+    assert np.allclose(val, 12.0225525)
 
     # close figure
     plt.close()
@@ -481,15 +481,15 @@ def test_specific_discharge_comprehensive():
     for col in ax.collections:
         if not isinstance(col, Quiver):
             raise AssertionError("Unexpected collection type")
-    assert np.sum(quiver.Umask) == 5
+    assert np.sum(quiver.Umask) == 7
     X = np.ma.masked_where(quiver.Umask, quiver.X)
     Y = np.ma.masked_where(quiver.Umask, quiver.Y)
     pos = X.sum() + Y.sum()
-    assert np.allclose(pos, -153.8341064453125)
+    assert np.allclose(pos, -152.0747652053833)
     U = np.ma.masked_where(quiver.Umask, quiver.U)
     V = np.ma.masked_where(quiver.Umask, quiver.V)
     val = U.sum() + V.sum()
-    assert np.allclose(val, -5.2501876516091235)
+    assert np.allclose(val, -3.3428158026088326)
 
     # close figure
     plt.close()

autotest/t076_modelgrid_thick.py

@@ -0,0 +1,205 @@
+"""
+Model grid thick method tests
+
+"""
+
+import numpy as np
+from flopy.discretization import StructuredGrid, VertexGrid, UnstructuredGrid
+
+
+def test_structured_thick():
+    nlay, nrow, ncol = 3, 2, 3
+    delc = 1.0 * np.ones(nrow, dtype=float)
+    delr = 1.0 * np.ones(ncol, dtype=float)
+    top = 10.0 * np.ones((nrow, ncol), dtype=float)
+    botm = np.zeros((nlay, nrow, ncol), dtype=float)
+    botm[0, :, :] = 5.0
+    botm[1, :, :] = 0.0
+    botm[2, :, :] = -5.0
+    grid = StructuredGrid(
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delc=delc,
+        delr=delr,
+        top=top,
+        botm=botm,
+    )
+    thick = grid.thick
+    assert np.allclose(thick, 5.0), "thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 10.0)
+    assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 5.0)
+    assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 2.5)
+    assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"
+
+    sat_thick = grid.saturated_thick(grid.botm)
+    assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
+
+    sat_thick = grid.saturated_thick(grid.botm - 100.0)
+    assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
+
+    return
+
+
+def test_vertices_thick():
+    nlay, ncpl = 3, 5
+    vertices = [
+        [0, 0.0, 3.0],
+        [1, 1.0, 3.0],
+        [2, 2.0, 3.0],
+        [3, 0.0, 2.0],
+        [4, 1.0, 2.0],
+        [5, 2.0, 2.0],
+        [6, 0.0, 1.0],
+        [7, 1.0, 1.0],
+        [8, 2.0, 1.0],
+        [9, 0.0, 0.0],
+        [10, 1.0, 0.0],
+    ]
+    iverts = [
+        [0, 0, 1, 4, 3],
+        [1, 1, 2, 5, 4],
+        [2, 3, 4, 7, 6],
+        [3, 4, 5, 8, 7],
+        [4, 6, 7, 10, 9],
+        [5, 0, 1, 4, 3],
+        [6, 1, 2, 5, 4],
+        [7, 3, 4, 7, 6],
+        [8, 4, 5, 8, 7],
+        [9, 6, 7, 10, 9],
+        [10, 0, 1, 4, 3],
+        [11, 1, 2, 5, 4],
+        [12, 3, 4, 7, 6],
+        [13, 4, 5, 8, 7],
+        [14, 6, 7, 10, 9],
+    ]
+    top = np.ones(ncpl, dtype=float) * 10.0
+    botm = np.zeros((nlay, ncpl), dtype=float)
+    botm[0, :] = 5.0
+    botm[1, :] = 0.0
+    botm[2, :] = -5.0
+    grid = VertexGrid(
+        nlay=nlay,
+        ncpl=ncpl,
+        vertices=vertices,
+        cell2d=iverts,
+        top=top,
+        botm=botm,
+    )
+    thick = grid.thick
+    assert np.allclose(thick, 5.0), "thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 10.0)
+    assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 5.0)
+    assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 2.5)
+    assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"
+
+    sat_thick = grid.saturated_thick(grid.botm)
+    assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
+
+    sat_thick = grid.saturated_thick(grid.botm - 100.0)
+    assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
+
+    return
+
+
+def test_unstructured_thick():
+    nlay = 3
+    ncpl = [5, 5, 5]
+    vertices = [
+        [0, 0.0, 3.0],
+        [1, 1.0, 3.0],
+        [2, 2.0, 3.0],
+        [3, 0.0, 2.0],
+        [4, 1.0, 2.0],
+        [5, 2.0, 2.0],
+        [6, 0.0, 1.0],
+        [7, 1.0, 1.0],
+        [8, 2.0, 1.0],
+        [9, 0.0, 0.0],
+        [10, 1.0, 0.0],
+    ]
+    iverts = [
+        [0, 0, 1, 4, 3],
+        [1, 1, 2, 5, 4],
+        [2, 3, 4, 7, 6],
+        [3, 4, 5, 8, 7],
+        [4, 6, 7, 10, 9],
+        [5, 0, 1, 4, 3],
+        [6, 1, 2, 5, 4],
+        [7, 3, 4, 7, 6],
+        [8, 4, 5, 8, 7],
+        [9, 6, 7, 10, 9],
+        [10, 0, 1, 4, 3],
+        [11, 1, 2, 5, 4],
+        [12, 3, 4, 7, 6],
+        [13, 4, 5, 8, 7],
+        [14, 6, 7, 10, 9],
+    ]
+    xcenters = [
+        0.5,
+        1.5,
+        0.5,
+        1.5,
+        0.5,
+    ]
+    ycenters = [
+        2.5,
+        2.5,
+        1.5,
+        1.5,
+        0.5,
+    ]
+    top = np.ones((nlay, 5), dtype=float)
+    top[0, :] = 10.0
+    top[1, :] = 5.0
+    top[2, :] = 0.0
+    botm = np.zeros((nlay, 5), dtype=float)
+    botm[0, :] = 5.0
+    botm[1, :] = 0.0
+    botm[2, :] = -5.0
+
+    grid = UnstructuredGrid(
+        vertices=vertices,
+        iverts=iverts,
+        xcenters=xcenters,
+        ycenters=ycenters,
+        ncpl=ncpl,
+        top=top.flatten(),
+        botm=botm.flatten(),
+    )
+
+    thick = grid.thick
+    assert np.allclose(thick, 5.0), "thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 10.0)
+    assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 5.0)
+    assert np.allclose(sat_thick, thick), "saturated thicknesses != 5."
+
+    sat_thick = grid.saturated_thick(grid.botm + 2.5)
+    assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"
+
+    sat_thick = grid.saturated_thick(grid.botm)
+    assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
+
+    sat_thick = grid.saturated_thick(grid.botm - 100.0)
+    assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
+
+    return
+
+
+if __name__ == "__main__":
+    test_unstructured_thick()
+    test_vertices_thick()
+    test_structured_thick()

autotest/t504_test.py

@@ -925,8 +925,11 @@ def test045_lake2tr():
         head_file = os.path.join(os.getcwd(), expected_head_file_a)
         head_new = os.path.join(run_folder, "lakeex2a.hds")
         assert pymake.compare_heads(
-            None, None, files1=head_file, files2=head_new,
-            htol=10.,
+            None,
+            None,
+            files1=head_file,
+            files2=head_new,
+            htol=10.0,
         )
 
     # change some settings
@@ -952,8 +955,11 @@ def test045_lake2tr():
         head_file = os.path.join(os.getcwd(), expected_head_file_b)
         head_new = os.path.join(save_folder, "lakeex2a.hds")
         assert pymake.compare_heads(
-            None, None, files1=head_file, files2=head_new,
-            htol=10.,
+            None,
+            None,
+            files1=head_file,
+            files2=head_new,
+            htol=10.0,
         )
 
 
@@ -1084,8 +1090,12 @@ def test027_timeseriestest():
         head_new = os.path.join(run_folder, "timeseriestest.hds")
         outfile = os.path.join(run_folder, "head_compare.dat")
         assert pymake.compare_heads(
-            None, None, files1=head_file, files2=head_new, outfile=outfile,
-            htol=10.,
+            None,
+            None,
+            files1=head_file,
+            files2=head_new,
+            outfile=outfile,
+            htol=10.0,
         )
 
     model = sim.get_model(model_name)
@@ -1108,8 +1118,11 @@ def test027_timeseriestest():
         head_file = os.path.join(os.getcwd(), expected_head_file_b)
         head_new = os.path.join(save_folder, "timeseriestest.hds")
         assert pymake.compare_heads(
-            None, None, files1=head_file, files2=head_new,
-            htol=10.,
+            None,
+            None,
+            files1=head_file,
+            files2=head_new,
+            htol=10.0,
         )
 
 

autotest/t505_test.py

@@ -2916,8 +2916,12 @@ def test028_sfr():
         head_new = os.path.join(run_folder, "test1tr.hds")
         outfile = os.path.join(run_folder, "head_compare.dat")
         assert pymake.compare_heads(
-            None, None, files1=head_file, files2=head_new, outfile=outfile,
-            htol=10.,
+            None,
+            None,
+            files1=head_file,
+            files2=head_new,
+            outfile=outfile,
+            htol=10.0,
         )
 
         # clean up