api: add mul interp mode

Author: mloubout

devito/finite_differences/derivative.py

@@ -583,7 +583,7 @@ def _eval_fd(self, expr, **kwargs):
         if self.method == 'RSFD':
             assert len(self.dims) == 1
             assert self.deriv_order[0] == 1
-            res = d45(expr, self.dims[0], x0=self.x0, expand=expand)
+            res = d45(expr, self.dims[0], x0=x0_deriv, expand=expand)
         elif len(self.dims) > 1:
             assert self.method == 'FD'
             res = cross_derivative(expr, self.dims, self.fd_order, self.deriv_order,
@@ -594,7 +594,7 @@ def _eval_fd(self, expr, **kwargs):
             res = generic_derivative(expr, self.dims[0], self.fd_order[0],
                                      self.deriv_order[0], weights=self.weights,
                                      side=self.side, matvec=self.transpose,
-                                     x0=self.x0, expand=expand)
+                                     x0=x0_deriv, expand=expand)
 
         # Step 4: Apply substitutions
         for e in self._ppsubs:

devito/finite_differences/differentiable.py

@@ -629,6 +629,51 @@ def _gather_for_diff(self):
 
         return self.func(*new_args, evaluate=False)
 
+    def _eval_at(self, func):
+        # No a basic a*b*c... expression, just defer to superclass
+        if any(isinstance(f, DifferentiableOp) for f in self.args):
+            return super()._eval_at(func)
+
+        # Split Derivative and Differentiable args
+        derivs, other = split(self.args, lambda e: isinstance(e, sympy.Derivative))
+
+        if derivs:
+            derivs = Differentiable._eval_at(self.func(*derivs), func)
+        else:
+            derivs = 1
+
+        if not other:
+            return derivs
+        elif len(other) > 1:
+            expr = self.func(*other)._gather_for_diff
+        else:
+            expr = other[0]
+
+        # Non differentiable expr (e.g., number)
+        if not isinstance(expr, Differentiable):
+            return self.func(derivs, expr)
+
+        # Build mapper for dimensions that need to be interpolated
+        mapper = {}
+        for d in self.dimensions:
+            try:
+                if self.indices_ref[d] is not func.indices_ref[d]:
+                    mapper[d] = func.indices_ref[d]
+            except KeyError:
+                pass
+
+        # Nothing to interpolate
+        if not mapper:
+            return super()._eval_at(func)
+
+        # Interpolate expr at the required indices
+        interp = expr.diff(*mapper.keys(), deriv_order=[0 for _ in mapper],
+                           fd_order=[self.interp_order for _ in mapper],
+                           x0=mapper)
+
+        # Return the full expression with Derivatives
+        return self.func(derivs, interp)
+
 
 class Pow(DifferentiableOp, sympy.Pow):
     _fd_priority = 0

devito/types/dense.py

@@ -1070,7 +1070,7 @@ def __fd_setup__(self):
 
     @cached_property
     def _fd_priority(self):
-        return 1 if self.staggered.on_node else 2
+        return 2.1 if self.staggered.on_node else 2
 
     def _eval_at(self, func):
         if self.staggered == func.staggered:
@@ -1491,7 +1491,7 @@ def __shape_setup__(cls, **kwargs):
 
     @cached_property
     def _fd_priority(self):
-        return 2.1 if self.staggered.on_node else 2.2
+        return 2.3 if self.staggered.on_node else 2.2
 
     @property
     def time_order(self):