From fd238c63e649b7463c7926bb17eed0eb34402ed4 Mon Sep 17 00:00:00 2001
From: Eric Barefoot <eabarefoot@gmail.com>
Date: Mon, 25 May 2020 20:12:11 -0500
Subject: [PATCH] fix: water and sediment routing bugs that caused inconsistent
 behavior

There were two issues with the refactoring. One was that water routing
developed a flat flow field because I wasn't adding each passing water
parcel to the Q field. The other was related to sediment parcel routing,
where I had transcribed the i and j of stepping incorrectly in the
jitted functions.
---
 pyDeltaRCM/sed_tools.py    |  6 +++---
 pyDeltaRCM/shared_tools.py | 17 +++++++++++------
 pyDeltaRCM/water_tools.py  |  2 +-
 3 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/pyDeltaRCM/sed_tools.py b/pyDeltaRCM/sed_tools.py
index 0c456c23..ece4ebc4 100644
--- a/pyDeltaRCM/sed_tools.py
+++ b/pyDeltaRCM/sed_tools.py
@@ -207,14 +207,14 @@ def sed_parcel(self, theta_sed, sed, px, py):
 
                 depoPart = self.Vp_res / 2 / self.dt / self.dx
 
-                px, py, self.qs[px, py] = shared_tools.partition_sand(self.qs[px, py], depoPart, py, px, dist, istep, jstep)
+                px, py, self.qs = shared_tools.partition_sand(self.qs, depoPart, py, px, dist, istep, jstep)
 
                 self.sand_dep_ero(px, py)
 
             if sed == 'mud':  # mud
 
-                px = px + istep
-                py = py + jstep
+                px = px + jstep
+                py = py + istep
 
                 self.mud_dep_ero(px, py)
 
diff --git a/pyDeltaRCM/shared_tools.py b/pyDeltaRCM/shared_tools.py
index 11ebed4a..e6d525e2 100644
--- a/pyDeltaRCM/shared_tools.py
+++ b/pyDeltaRCM/shared_tools.py
@@ -62,13 +62,14 @@ def get_filtered_weight(weight, px, depth_ind, cell_type_ind, dry_depth):
 def partition_sand(qs, depoPart, py, px, dist, istep, jstep):
     if dist > 0:
         # deposition in current cell
-        qs = qs + depoPart
-    px = px + istep
-    py = py + jstep
+        qs[px, py] += depoPart
+
+    px = px + jstep
+    py = py + istep
 
     if dist > 0:
         # deposition in downstream cell
-        qs = qs + depoPart
+        qs[px, py] += depoPart
     return px, py, qs
 
 
@@ -86,13 +87,17 @@ def get_steps(new_cells, iwalk, jwalk):
 
 @njit
 def update_dirQfield(qfield, dist, inds, astep, dirstep):
-    qfield[inds[astep]] += dirstep[astep] / dist[astep]
+    for i, ii in enumerate(inds):
+        if astep[i]:
+            qfield[ii] += dirstep[i] / dist[i]
     return qfield
 
 
 @njit
 def update_absQfield(qfield, dist, inds, astep, Qp_water, dx):
-    qfield[inds[astep]] += Qp_water / dx / 2
+    for i, ii in enumerate(inds):
+        if astep[i]:
+            qfield[ii] += Qp_water / dx / 2
     return qfield
 
 
diff --git a/pyDeltaRCM/water_tools.py b/pyDeltaRCM/water_tools.py
index 41579cde..9ea5456a 100644
--- a/pyDeltaRCM/water_tools.py
+++ b/pyDeltaRCM/water_tools.py
@@ -89,7 +89,7 @@ def run_water_iteration(self):
 
             current_inds, self.looped, self.free_surf_flag = shared_tools.check_for_loops(
                 self.indices,
-                current_inds,
+                next_index,
                 iter,
                 self.L0,
                 self.looped,