Repository: devitocodes/devito
Issue: Error with wrong indices in loops

The following MFE presents some cases in which the indexes used to access functions are not defined correctly. All these cases produce more than one time-iteration loop and that's when the error appears. 

```python
from devito import (Eq, Operator, Grid, TimeFunction)

grid = Grid(shape=(10,10))
x, y = grid.dimensions
t = grid.stepping_dim

f = TimeFunction(name="f", grid=grid, space_order=2, time_order=2) 
g = TimeFunction(name="g", grid=grid, space_order=2, time_order=2)

eqn = []
eqn.append( Eq( f.forward, f + 1 ) )

# Case 1 - Compilation error (t0 undeclared on the second loop, t0+2 should be replaced by t2)
# eqn.append( Eq( g.forward, f[t+2,x,y] ) )

# Case 2 - Segfault (t0 is declared, but t0+3 attempts to access wrong index, t0+3 should be replaced by t2)
eqn.append( Eq( g, f[t+3,x,y] ) ) 

# Case 3 - No error (t0+2 should be replaced by t2)
# eqn.append( Eq( g, f[t+2,x,y] ) )

op = Operator(eqn)

op.apply(time_M=10)
```

It seems substitutions of `SteppingDimensions` by `ModuloDimensions` in  `_lower_stepping_dims(iet)` (devito/ir/iet/scheduler.py) are being done in a misleading manner: All expressions in `iet` have their indices replaced with `uindices` (eg. `t+2` -> `t0+2`) available in the first iteration tree, which in these examples is the first time-iteration loop. That way, when we get to the second time-iteration loop, indices that were supposed to be replaced (eg. `t+2` -> `t2`) are not found because they have already been replaced (eg. `t+2` -> `t0+2`) and so they remain that way in the code, resulting in compilation or execution errors. This is because the second time-iteration loop may not even have that `t0` indice (Case 1) or if it has `t0` it can attempt to access wrong index (Case 2) since indexes are miscalculated (eg. `t0+2 != t2`, because `t0 = (time)%(3)`, `t2 = (time+2)%(3)`, so `t0+2` could reach 4, out-of-bounds, as soon as t0==2). Also, that could even not produce any errors at all, but still mess up the numerical result (Case 3).

Changed files:
devito/ir/iet/scheduler.py, tests/test_operator.py

diff --git a/devito/ir/iet/scheduler.py b/devito/ir/iet/scheduler.py
index 693d6adf6..03f37ea62 100644
--- a/devito/ir/iet/scheduler.py
+++ b/devito/ir/iet/scheduler.py
@@ -2,7 +2,7 @@ from collections import OrderedDict
 
 from devito.ir.iet import (Expression, Increment, Iteration, List, Conditional,
                            Section, HaloSpot, ExpressionBundle, FindNodes, FindSymbols,
-                           XSubs)
+                           XSubs, Transformer)
 from devito.symbolics import IntDiv, xreplace_indices
 from devito.tools import as_mapper, timed_pass
 from devito.types import ConditionalDimension
@@ -69,6 +69,7 @@ def _lower_stepping_dims(iet):
 
     u[t1, x] = u[t0, x] + 1
     """
+    subs = {}
     for i in FindNodes(Iteration).visit(iet):
         if not i.uindices:
             # Be quick: avoid uselessy reconstructing nodes
@@ -80,13 +81,15 @@ def _lower_stepping_dims(iet):
         # two different calls to `xreplace`
         mindices = [d for d in i.uindices if d.is_Modulo]
         groups = as_mapper(mindices, lambda d: d.modulo)
+        root = i
         for k, v in groups.items():
             mapper = {d.origin: d for d in v}
-            rule = lambda i: i.function.is_TimeFunction and i.function._time_size == k
-            replacer = lambda i: xreplace_indices(i, mapper, rule)
-            iet = XSubs(replacer=replacer).visit(iet)
+            rule = lambda e: e.function.is_TimeFunction and e.function._time_size == k
+            replacer = lambda e: xreplace_indices(e, mapper, rule)
+            root = XSubs(replacer=replacer).visit(root)
+        subs[i] = root
 
-    return iet
+    return Transformer(subs).visit(iet)
 
 
 def _lower_conditional_dims(iet):
diff --git a/tests/test_operator.py b/tests/test_operator.py
index 563044cba..70ed19c4f 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -213,6 +213,42 @@ class TestCodeGen(object):
         assert(i.indices[i.function._time_position].modulo == exp_mods[i.function.name]
                for i in flatten(retrieve_indexed(i) for i in exprs))
 
+    def test_lower_stepping_dims_with_mutiple_iterations(self):
+        """
+        Test lowering SteppingDimensions for a time dimension with
+        more than one iteration loop with different ModuloDimensions.
+        MFE for issue #1486
+        """
+        grid = Grid(shape=(4, 4))
+
+        f = Function(name="f", grid=grid, space_order=4)
+        g = Function(name="g", grid=grid, space_order=4)
+        h = TimeFunction(name="h", grid=grid, space_order=4, time_order=2)
+
+        f.data[:] = 0.0
+        h.data[:] = 0.0
+
+        eqn = [Eq(f, h + 1), Eq(g, f),
+               Eq(h.forward, h + g + 1)]
+
+        op = Operator(eqn)
+
+        for iter in [i for i in FindNodes(Iteration).visit(op) if i.dim.is_Time]:
+            exprtimeindices = set([a.indices[a.function._time_position] for
+                                   expr in FindNodes(Expression).visit(iter) for
+                                   a in retrieve_indexed(expr.expr) if
+                                   isinstance(a.function, TimeFunction)])
+            # Check if iteration time indices match with expressions time indices
+            assert (exprtimeindices == set(iter.uindices))
+            # Check if expressions time indices are modulo dimensions
+            assert(all([i.is_Modulo for i in exprtimeindices]))
+
+        op.apply(time_M=10)
+
+        assert np.all(h.data[0, :] == 18)
+        assert np.all(h.data[1, :] == 20)
+        assert np.all(h.data[2, :] == 22)
+
     @skipif('device')
     def test_timedlist_wraps_time_if_parallel(self):
         """