closer

modflowpy · Sep 29, 2023 · ec9d124 · ec9d124
1 parent 68b83d6
commit ec9d124
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 13 deletions.
diff --git a/autotest/test_postprocessing.py b/autotest/test_postprocessing.py
@@ -39,7 +39,7 @@ def mf6_freyberg_path(example_data_path):
 
 @pytest.mark.mf6
 @requires_exe("mf6")
-def test_get_structured_faceflows_2d_right_lower(function_tmpdir):
+def test_get_structured_faceflows_issue_1911(function_tmpdir):
     """
     Reproduce https://github.com/modflowpy/flopy/issues/1911
     """
@@ -211,9 +211,15 @@ def test_get_structured_faceflows_2d_right_lower(function_tmpdir):
 @pytest.mark.parametrize(
     "nlay, nrow, ncol",
     [
+        # extended in 1 dimension
         [5, 1, 1],
         [1, 5, 1],
         [1, 1, 5],
+        # 2D
+        [5, 5, 1],
+        [1, 5, 5],
+        [5, 1, 5],
+        # 3D
     ],
 )
 @pytest.mark.mf6
@@ -309,9 +315,10 @@ def test_get_structured_faceflows_1d(function_tmpdir, nlay, nrow, ncol):
         verbose=True,
     )
 
-    assert np.any(frf) == bool(ncol > 1)
-    assert np.any(fff) == bool(nrow > 1)
-    assert np.any(flf) == bool(nlay > 1)
+    # expect nonzero flows only in extended (>1 cell) dimensions
+    assert np.any(frf) == (ncol > 1)
+    assert np.any(fff) == (nrow > 1)
+    assert np.any(flf) == (nlay > 1)
 
 
 @pytest.mark.skip(reason="todo")
@@ -369,15 +376,15 @@ def test_get_structured_faceflows_freyberg(
 
     # get freyberg mf2005 output
     mf2005_cbc = flopy.utils.CellBudgetFile(mf2005_ws / "freyberg.cbc")
-    # mf2005_spdis = mf2005_cbc.get_data(text="DATA-SPDIS")[0]
     mf2005_frf, mf2005_fff = (
         mf2005_cbc.get_data(text="FLOW RIGHT FACE", full3D=True)[0],
         mf2005_cbc.get_data(text="FLOW FRONT FACE", full3D=True)[0],
     )
 
-    assert mf2005_frf.shape == mf2005_fff.shape == mf6_head.shape
-    assert np.allclose(mf6_frf, mf2005_frf)
-    assert np.allclose(mf6_fff, mf2005_fff)
+    # todo compare mf2005 faceflows with converted mf6 faceflows (debug disagreement)
+    # assert mf2005_frf.shape == mf2005_fff.shape == mf6_head.shape
+    # assert np.allclose(mf6_frf, mf2005_frf)
+    # assert np.allclose(mf6_fff, mf2005_fff)
 
     Qx, Qy, Qz = get_specific_discharge(
         (mf6_frf, mf6_fff, mf6_flf),

diff --git a/flopy/mf6/utils/postprocessing.py b/flopy/mf6/utils/postprocessing.py
@@ -10,8 +10,9 @@ def get_face(m, n, nlay, nrow, ncol):
 
     Notes
     -----
-    For visual intuition see
+    For visual intuition in 2 dimensions
     https://stackoverflow.com/a/16330162/6514033
+    is helpful
 
     Parameters
     ----------
@@ -31,8 +32,15 @@ def get_face(m, n, nlay, nrow, ncol):
     face : int
         0: right, 1: front, 2: lower
     """
-    if m - 1 == n:
-        return 0
+
+    if m - n == 1:
+        # handle 1D cases
+        if nrow == 1 and ncol == 1:
+            return 2
+        elif nlay == 1 and ncol == 1:
+            return 1
+        else:
+            return 0
     elif m - n == nrow * ncol:
         return 2
     else:
@@ -47,7 +55,7 @@ def get_structured_faceflows(
     nlay=None,
     nrow=None,
     ncol=None,
-    verbose=False
+    verbose=False,
 ):
     """
     Get the face flows for the flow right face, flow front face, and
@@ -94,7 +102,13 @@ def get_structured_faceflows(
         ia, ja = grb.ia, grb.ja
         nlay, nrow, ncol = grb.nlay, grb.nrow, grb.ncol
     else:
-        if ia is None or ja is None or nlay is None or nrow is None or ncol is None:
+        if (
+            ia is None
+            or ja is None
+            or nlay is None
+            or nrow is None
+            or ncol is None
+        ):
             raise ValueError(
                 "ia, ja, nlay, nrow, and ncol must be specified if the MODFLOW 6"
                 "binary grid file name is not specified."