diff --git a/devito/ir/equations/algorithms.py b/devito/ir/equations/algorithms.py
index 2606d2dd2dc..6911082dff3 100644
--- a/devito/ir/equations/algorithms.py
+++ b/devito/ir/equations/algorithms.py
@@ -39,7 +39,13 @@ def handle_indexed(indexed):
         # passes, so they can be ignored here
         relation = tuple(d for d in relation if not d.is_Stencil)
 
-        return relation
+        # Make sure indices appearance order satisfies "parent" order. Foe example
+        # case where f[rx, ps] where ps is a parent of `rx` would lead to ordering issues
+        # as the order of relation is assumed to be mandatory
+        rel = {(d.root, d) for d in relation if d.is_Derived and d.root in relation}
+        relation = PartialOrderTuple(relation, relations=rel)
+
+        return tuple(relation)
 
     if isinstance(expr.implicit_dims, IgnoreDimSort):
         relations = set()
@@ -74,7 +80,7 @@ def handle_indexed(indexed):
     # obtain the desired ordering `(time, x, xi)`. W/o `(time, x)`, the ordering
     # `(x, time, xi)` might be returned instead, which would be non-sense
     for i in relations:
-        dims = [di for d in i for di in (d.index, d)]
+        dims = [di for d in i for di in (d.root, d)]
         implicit_relations.update({tuple(filter_ordered(dims))})
 
     ordering = PartialOrderTuple(extra, relations=(relations | implicit_relations))
diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py
index e082848ffca..3f1ad5e3b63 100644
--- a/devito/operations/interpolators.py
+++ b/devito/operations/interpolators.py
@@ -176,11 +176,7 @@ def _interp_idx(self, variables, implicit_dims=None):
         """
         mapper = {}
         pos = self.sfunction._position_map.values()
-        # Temporaries for the position
-        temps = self._positions(implicit_dims)
 
-        # Coefficient symbol expression
-        temps.extend(self._coeff_temps(implicit_dims))
         for ((di, d), rd, p) in zip(enumerate(self._gdims), self._rdim, pos):
             # Add conditional to avoid OOB
             lb = sympy.And(rd + p >= d.symbolic_min - self.r, evaluate=False)
@@ -188,10 +184,17 @@ def _interp_idx(self, variables, implicit_dims=None):
             cond = sympy.And(lb, ub, evaluate=False)
             mapper[d] = ConditionalDimension(rd.name, rd, condition=cond, indirect=True)
 
+        # Temporaries for the position
+        temps = self._positions(implicit_dims)
+
+        # Coefficient symbol expression
+        temps.extend(self._coeff_temps(implicit_dims))
+
         # Substitution mapper for variables
         idx_subs = {v: v.subs({k: c - v.origin.get(k, 0) + p
                     for ((k, c), p) in zip(mapper.items(), pos)})
                     for v in variables}
+        idx_subs.update(dict(zip(self._rdim, mapper.values())))
 
         return idx_subs, temps