fix(lgrutil): corrected bug in child to parent vertical exchanges (#948)

autotest/t063_test_lgrutil.py

@@ -45,28 +45,37 @@ def test_lgrutil():
     assert ac[0, 0] == 6
     assert ac[-1, -1] == 18
 
-    # top/bottom
-    top, botm = lgr.get_top_botm()
-    assert top.shape == (9, 9)
-    assert botm.shape == (2, 9, 9)
-    assert top.min() == top.max() == 100.
-    assert botm.min() == botm.max() == 0.
+    # child top/bottom
+    topc, botmc = lgr.get_top_botm()
+    assert topc.shape == (9, 9)
+    assert botmc.shape == (2, 9, 9)
+    assert topc.min() == topc.max() == 100.
+    errmsg = '{} /= {}'.format(botmc[:, 0, 0], np.array(botmp[:2]))
+    assert np.allclose(botmc[:, 0, 0], np.array(botmp[:2])), errmsg
 
     # exchange data
     exchange_data = lgr.get_exchange_data(angldegx=True, cdist=True)
+
+
     ans1 = [(0, 1, 0), (0, 0, 0), 1, 50.0, 16.666666666666668,
             33.333333333333336, 0.0, 354.33819375782156]
-    ans2 = [(1, 4, 3), (1, 8, 8), 1, 50.0, 16.666666666666668,
-            33.333333333333336, 90.0, 354.3381937578216]
-    assert exchange_data[0] == ans1
-    assert exchange_data[-1] == ans2
-    assert len(exchange_data) == 72
+    errmsg = '{} /= {}'.format(ans1, exchange_data[0])
+    assert exchange_data[0] == ans1, errmsg
+
+    ans2 = [(2, 3, 3), (1, 8, 8), 0, 50.0, 50,
+            1111.1111111111113, 180., 100.]
+    errmsg = '{} /= {}'.format(ans2, exchange_data[-1])
+    assert exchange_data[-1] == ans2, errmsg
+
+    errmsg = 'exchanges should be 71 horizontal plus 81 vertical'
+    assert len(exchange_data) == 72 + 81, errmsg
 
     # list of parent cells connected to a child cell
     assert lgr.get_parent_connections(0, 0, 0) == [((0, 1, 0), -1),
                                                    ((0, 0, 1), 2)]
     assert lgr.get_parent_connections(1, 8, 8) == [((1, 3, 4), 1),
-                                                   ((1, 4, 3), -2)]
+                                                   ((1, 4, 3), -2),
+                                                   ((2, 3, 3), -3)]
 
     return
 

flopy/utils/lgrutil.py

@@ -66,6 +66,14 @@ def __init__(self, nlayp, nrowp, ncolp, delrp, delcp, topp, botmp,
         self.ncpp = ncpp
         self.ncppl = Util2d(m, (nlayp,), np.int, ncppl, 'ncppl').array
 
+        # calculate ibcl which is the bottom child layer (one based) in each
+        # parent layer
+        self.ibcl = np.zeros(self.nlayp, dtype=np.int)
+        self.ibcl[0] = self.ncppl[0]
+        for k in range(1, self.nlayp):
+            if self.ncppl[k] > 0:
+                self.ibcl[k] = self.ibcl[k - 1] + self.ncppl[k]
+
         # parent lower left
         self.xllp = xllp
         self.yllp = yllp
@@ -152,8 +160,8 @@ def get_top_botm(self):
                 for kp in range(self.nplbeg, self.nplend + 1):
                     top = pbotm[kp, ip, jp]
                     bot = pbotm[kp + 1, ip, jp]
-                    dz = (top - bot) / self.ncppl[kp - 1]
-                    for _ in range(self.ncppl[kp - 1]):
+                    dz = (top - bot) / self.ncppl[kp]
+                    for _ in range(self.ncppl[kp]):
                         botm[kc, icrowstart:icrowend,
                         iccolstart: iccolend] = botm[kc - 1,
                                                 icrowstart:icrowend,
@@ -273,7 +281,7 @@ def get_parent_connections(self, kc, ic, jc):
         # parent cell to top is not possible
 
         # parent cell to bottom
-        if kc + 1 == self.ncppl[kp]:
+        if kc + 1 == self.ibcl[kp]:
             if kp + 1 < self.nlayp:
                 if self.idomain[kp + 1, ip, jp] != 0:
                     parentlist.append(((kp + 1, ip, jp), -3))
@@ -335,6 +343,12 @@ def get_exchange_data(self, angldegx=False, cdist=False):
                             continue
 
                         # horizontal or vertical connection
+                        # 1 if a child cell horizontally connected to a parent
+                        #   cell
+                        # 2 if more than one child cells horizontally connected
+                        #   to parent cell
+                        # 0 if a vertical connection
+
                         ihc = 1
                         if self.ncppl[kp] > 1:
                             ihc = 2