From 0616db5f158a97d1b38544d993d866c19a672a80 Mon Sep 17 00:00:00 2001
From: langevin-usgs <langevin@usgs.gov>
Date: Fri, 21 Jan 2022 17:54:05 -0600
Subject: [PATCH] fix(postprocessing): get_structured_faceflows fix to support
 3d models (#1333)

* function was fixed to support multi-layer models
* test added to ensure face flow calculations work
---
 autotest/t029_test.py             | 40 +++++++++++++++++++++++++++++--
 flopy/mf6/utils/postprocessing.py |  2 +-
 2 files changed, 39 insertions(+), 3 deletions(-)

diff --git a/autotest/t029_test.py b/autotest/t029_test.py
index e1128412df..1e65c23236 100644
--- a/autotest/t029_test.py
+++ b/autotest/t029_test.py
@@ -326,6 +326,41 @@ def test_flowja_residuals():
     return
 
 
+def test_structured_faceflows_3d():
+    model_ws = f"{base_dir}_test_faceflows_3d"
+    test_setup = FlopyTestSetup(verbose=True, test_dirs=model_ws)
+    name = 'mymodel'
+    sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=model_ws, exe_name='mf6')
+    tdis = flopy.mf6.ModflowTdis(sim)
+    ims = flopy.mf6.ModflowIms(sim)
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
+    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=3, nrow=10, ncol=10, top=0, botm=[-1, -2, -3])
+    ic = flopy.mf6.ModflowGwfic(gwf)
+    npf = flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True)
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 0, 0), 1.],
+                                                           [(0, 9, 9), 0.]])
+    budget_file = name + '.bud'
+    head_file = name + '.hds'
+    oc = flopy.mf6.ModflowGwfoc(gwf,
+                                budget_filerecord=budget_file,
+                                head_filerecord=head_file,
+                                saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')])
+    sim.write_simulation()
+    sim.run_simulation()
+
+    head = gwf.output.head().get_data()
+    bud = gwf.output.budget()
+    flowja = bud.get_data(text="FLOW-JA-FACE")[0]
+    frf, fff, flf = flopy.mf6.utils.get_structured_faceflows(
+        flowja,
+        grb_file=os.path.join(model_ws, "mymodel.dis.grb"),
+    )
+    assert frf.shape == head.shape, f"frf.shape {frf.shape} != head.shape {head.shape}"
+    assert fff.shape == head.shape, f"frf.shape {frf.shape} != head.shape {head.shape}"
+    assert flf.shape == head.shape, f"frf.shape {frf.shape} != head.shape {head.shape}"
+    return
+
+
 def test_faceflows_empty():
     flowja = np.zeros(10, dtype=np.float64)
     with pytest.raises(ValueError):
@@ -375,5 +410,6 @@ def test_residuals_iaempty():
     # test_mfgrddisv_modelgrid()
     # test_mfgrddisu_MfGrdFile()
     # test_mfgrddisu_modelgrid()
-    test_faceflows()
-    test_flowja_residuals()
+    # test_faceflows()
+    test_structured_faceflows_3d()
+    # test_flowja_residuals()
diff --git a/flopy/mf6/utils/postprocessing.py b/flopy/mf6/utils/postprocessing.py
index bf4bee94bc..2122f79d66 100644
--- a/flopy/mf6/utils/postprocessing.py
+++ b/flopy/mf6/utils/postprocessing.py
@@ -60,7 +60,7 @@ def get_structured_faceflows(
     shape = (grb.nlay, grb.nrow, grb.ncol)
     frf = np.zeros(shape, dtype=float).flatten()
     fff = np.zeros(shape, dtype=float).flatten()
-    flf = np.zeros(shape, dtype=float)
+    flf = np.zeros(shape, dtype=float).flatten()
 
     # fill flow terms
     vmult = [-1.0, -1.0, -1.0]