dsl: Handle nested derivatives that have specified fd_order

Author: JDBetteridge

devito/finite_differences/derivative.py

@@ -149,16 +149,31 @@ def __new__(cls, expr, *dims, **kwargs):
             else:
                 dcounter[d] += o
 
-        # Default finite difference orders depending on input dimension (.dt or .dx)
         # It's possible that the expr is a `sympy.Number` at this point, which
         # has derivative 0, unless we're taking a 0th derivative.
         if isinstance(expr, sympy.Number):
             if any(dcounter.values()):
                 return 0
             else:
                 return expr
+
+        # Use `fd_order` if specified
+        fd_order = kwargs.get('fd_order')
+        if fd_order is not None:
+            _fd_order_specified = True
+            # If `fd_order` is specified collect these together
+            fcounter = defaultdict(int)
+            for d, o in zip(dims, as_tuple(fd_order)):
+                fcounter[d] += o
+            for d, o in dcounter.items():
+                order = expr.time_order if getattr(d, 'is_Time', False) else expr.space_order
+                if o > order:
+                    raise ValueError(f'Function does not support {d} derivative of order {o}')
+            fd_order = fcounter.values()
         else:
-            default_fdo = tuple([
+            _fd_order_specified = False
+            # Default finite difference orders depending on input dimension (.dt or .dx)
+            fd_order = tuple([
                 expr.time_order
                 if getattr(d, 'is_Time', False)
                 else expr.space_order
@@ -174,8 +189,9 @@ def __new__(cls, expr, *dims, **kwargs):
         obj = Differentiable.__new__(cls, expr, *derivatives)
         obj._dims = tuple(dcounter.keys())
 
+        obj._fd_order_specified = _fd_order_specified
         obj._fd_order = DimensionTuple(
-            *as_tuple(kwargs.get('fd_order', default_fdo)),
+            *as_tuple(fd_order),
             getters=obj._dims
         )
         obj._deriv_order = DimensionTuple(
@@ -533,7 +549,11 @@ def _eval_expand_nest(self, **hints):
             # all kwargs from the self object
             # TODO: This dictionary merge needs to be a lot better
             # EG: Don't actually expand if derivatives are incompatible
-            new_kwargs = {'deriv_order': tuple(chain(self.deriv_order, expr.deriv_order))}
+            new_kwargs = {
+                'deriv_order': tuple(chain(self.deriv_order, expr.deriv_order))
+            }
+            if self._fd_order_specified or expr._fd_order_specified:
+                new_kwargs.update({'fd_order': tuple(chain(self.fd_order, expr.fd_order))})
             return self.func(new_expr, *new_dims, **new_kwargs)
         else:
             return self
@@ -564,3 +584,16 @@ def _eval_expand_add(self, **hints):
                 return self.func(dep, *self.args[1:])
         else:
             return self
+
+    def _eval_expand_product_rule(self, **hints):
+        ''' Expands products, of functions of the dependent variable
+        `Derivative(f(x)·g(x), x)
+            --> Derivative(f(x), x)·g(x) + f(x)·Derivative(g(x), x)`
+        This is only implemented for first derivatives with an arbitrary number
+        of multiplicands and second derivatives with two multiplicands. The
+        resultant expression for higher derivatives and mixed derivatives is much
+        more difficult to implement.
+        '''
+        # expr = self.args[0]
+        # if isinstance(expr, sympy.Mul):
+        raise NotImplementedError('Product rule expansion has not been written')

tests/test_derivatives.py

@@ -1051,19 +1051,19 @@ def test_deriv_order(self):
         assert du11 == du2
         assert du11.deriv_order == du2.deriv_order
 
-    @pytest.mark.xfail(raises=ValueError)
     def test_wrong_deriv_order(self):
         ''' Check an exception is raises with incompatible arguments
         '''
-        _ = Derivative(self.u, self.x, deriv_order=(2, 4))
+        with pytest.raises(ValueError):
+            _ = Derivative(self.u, self.x, deriv_order=(2, 4))
 
-    @pytest.mark.xfail(raises=ValueError)
     def test_no_derivative(self):
-        _ = Derivative(sympify(-1))
+        with pytest.raises(ValueError):
+            _ = Derivative(sympify(-1))
 
-    @pytest.mark.xfail(raises=ValueError)
     def test_no_dimension(self):
-        _ = Derivative(sympify(-1), deriv_order=0)
+        with pytest.raises(ValueError):
+            _ = Derivative(sympify(-1), deriv_order=0)
 
     def test_constant(self):
         ''' Check constant derivative is zero for non-0th order derivatives
@@ -1144,3 +1144,41 @@ def test_expand_nest(self):
         '''
         expanded = 4*Derivative(self.u, self.x) - 5*Derivative(self.u, (self.x, 2))
         assert self.b.expand(add=True, nest=True) == expanded
+
+    def test_nested_orders(self):
+        ''' Check nested expansion results in correct derivative and fd order
+        '''
+        # Default fd_order
+        du22 = Derivative(Derivative(self.u, (self.x, 2)), (self.x, 2))
+        du22_expanded = du22.expand(nest=True)
+        du4 = Derivative(self.u, (self.x, 4))
+        assert du22_expanded == du4
+        assert du22_expanded.deriv_order == du4.deriv_order
+        assert du22_expanded.fd_order == du4.fd_order
+
+        # Specified fd_order
+        du22 = Derivative(
+            Derivative(self.u, (self.x, 2), fd_order=2),
+            (self.x, 2),
+            fd_order=2
+        )
+        du22_expanded = du22.expand(nest=True)
+        du4 = Derivative(self.u, (self.x, 4), fd_order=4)
+        assert du22_expanded == du4
+        assert du22_expanded.deriv_order == du4.deriv_order
+        assert du22_expanded.fd_order == du4.fd_order
+
+        # Specified fd_order greater than the space order
+        ## When no order specified
+        du44 = Derivative(Derivative(self.u, (self.x, 4)), (self.x, 4))
+        with pytest.raises(ValueError):
+            du44_expanded = du44.expand(nest=True)
+
+        ## When order specified is too large
+        du44 = Derivative(
+            Derivative(self.u, (self.x, 4), fd_order=4),
+            (self.x, 4),
+            fd_order=4
+        )
+        with pytest.raises(ValueError):
+            du44_expanded = du44.expand(nest=True)