api: Add shift and fd order option to all FD operators:

.github/workflows/tutorials.yml

@@ -30,7 +30,6 @@ jobs:
       matrix:
         name: [
           tutos-ubuntu-gcc-py39,
-          tutos-osx-gcc-py39,
           tutos-osx-clang-py39,
           tutos-docker-gcc-py39
           ]
@@ -41,11 +40,6 @@ jobs:
             compiler: gcc
             language: "openmp"
 
-          - name: tutos-osx-gcc-py39
-            os: macos-latest
-            compiler: gcc-11
-            language: "openmp"
-
           - name: tutos-osx-clang-py39
             os: macos-latest
             compiler: clang

devito/finite_differences/differentiable.py

@@ -269,21 +269,76 @@ def laplace(self):
         Generates a symbolic expression for the Laplacian, the second
         derivative w.r.t all spatial Dimensions.
         """
+        return self.laplacian()
+
+    def laplacian(self, shift=None, order=None):
+        """
+        Laplacian of the Differentiable with shifted derivatives and custom
+        FD order.
+
+        Each second derivative is left-right (i.e D^T D with D the first derivative ):
+        `(self.dx(x0=dim+shift*dim.spacing,
+                  fd_order=order)).dx(x0=dim-shift*dim.spacing, fd_order=order)`
+
+        Parameters
+        ----------
+        shift: Number, optional, default=None
+            Shift for the center point of the derivative in number of gridpoints
+        order: int, optional, default=None
+            Discretization order for the finite differences.
+            Uses `func.space_order` when not specified
+        """
+        order = order or self.space_order
         space_dims = [d for d in self.dimensions if d.is_Space]
+        shift_x0 = make_shift_x0(shift, (len(space_dims),))
         derivs = tuple('d%s2' % d.name for d in space_dims)
-        return Add(*[getattr(self, d) for d in derivs])
+        return Add(*[getattr(self, d)(x0=shift_x0(shift, space_dims[i], None, i),
+                                      fd_order=order)
+                     for i, d in enumerate(derivs)])
+
+    def div(self, shift=None, order=None):
+        """
+        Divergence of the input Function.
+
+        Parameters
+        ----------
+        func : Function or TensorFunction
+            Function to take the divergence of
+        shift: Number, optional, default=None
+            Shift for the center point of the derivative in number of gridpoints
+        order: int, optional, default=None
+            Discretization order for the finite differences.
+            Uses `func.space_order` when not specified
 
-    def div(self, shift=None):
+        """
         space_dims = [d for d in self.dimensions if d.is_Space]
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
-        return Add(*[getattr(self, 'd%s' % d.name)(x0=shift_x0(shift, d, None, i))
+        order = order or self.space_order
+        return Add(*[getattr(self, 'd%s' % d.name)(x0=shift_x0(shift, d, None, i),
+                                                   fd_order=order)
                      for i, d in enumerate(space_dims)])
 
-    def grad(self, shift=None):
+    def grad(self, shift=None, order=None):
+        """
+        Gradient of the input Function.
+
+        Parameters
+        ----------
+        func : Function or TensorFunction
+            Function to take the gradient of
+        shift: Number, optional, default=None
+            Shift for the center point of the derivative in number of gridpoints
+        order: int, optional, default=None
+            Discretization order for the finite differences.
+            Uses `func.space_order` when not specified
+
+        """
         from devito.types.tensor import VectorFunction, VectorTimeFunction
         space_dims = [d for d in self.dimensions if d.is_Space]
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
-        comps = [getattr(self, 'd%s' % d.name)(x0=shift_x0(shift, d, None, i))
+        order = order or self.space_order
+        comps = [getattr(self, 'd%s' % d.name)(x0=shift_x0(shift, d, None, i),
+                                               fd_order=order)
                  for i, d in enumerate(space_dims)]
         vec_func = VectorTimeFunction if self.is_TimeDependent else VectorFunction
         return vec_func(name='grad_%s' % self.name, time_order=self.time_order,

devito/finite_differences/operators.py

@@ -1,55 +1,81 @@
-def div(func, shift=None):
+def div(func, shift=None, order=None):
     """
     Divergence of the input Function.
 
     Parameters
     ----------
     func : Function or TensorFunction
+        Function to take the divergence of
+    shift: Number, optional, default=None
+        Shift for the center point of the derivative in number of gridpoints
+    order: int, optional, default=None
+        Discretization order for the finite differences.
+        Uses `func.space_order` when not specified
+
     """
     try:
-        return func.div(shift=shift)
+        return func.div(shift=shift, order=order)
     except AttributeError:
         return 0
 
 
-def grad(func, shift=None):
+def grad(func, shift=None, order=None):
     """
     Gradient of the input Function.
 
     Parameters
     ----------
-    func : Function or VectorFunction
+    func : Function or TensorFunction
+        Function to take the gradient of
+    shift: Number, optional, default=None
+        Shift for the center point of the derivative in number of gridpoints
+    order: int, optional, default=None
+        Discretization order for the finite differences.
+        Uses `func.space_order` when not specified
+
     """
     try:
-        return func.grad(shift=shift)
+        return func.grad(shift=shift, order=order)
     except AttributeError:
         raise AttributeError("Gradient not supported for class %s" % func.__class__)
 
 
-def curl(func):
+def curl(func, shift=None, order=None):
     """
-    Curl of the input Function.
+    Curl of the input Function. Only supported for VectorFunction
 
     Parameters
     ----------
     func : VectorFunction
+        VectorFunction to take curl of
+    shift: Number, optional, default=None
+        Shift for the center point of the derivative in number of gridpoints
+    order: int, optional, default=None
+        Discretization order for the finite differences.
+        Uses `func.space_order` when not specified
     """
     try:
-        return func.curl
+        return func.curl(shift=shift, order=order)
     except AttributeError:
         raise AttributeError("Curl only supported for 3D VectorFunction")
 
 
-def laplace(func):
+def laplace(func, shift=None, order=None):
     """
     Laplacian of the input Function.
 
     Parameters
     ----------
-    func : Function or TensorFunction
+    func : VectorFunction
+        VectorFunction to take laplacian of
+    shift: Number, optional, default=None
+        Shift for the center point of the derivative in number of gridpoints
+    order: int, optional, default=None
+        Discretization order for the finite differences.
+        Uses `func.space_order` when not specified
     """
     try:
-        return func.laplace
+        return func.laplacian(shift=shift, order=order)
     except AttributeError:
         return 0
 
@@ -61,6 +87,10 @@ def diag(func, size=None):
     Parameters
     ----------
     func : Differentiable or scalar
+        Symbolic object to set the diagonal to
+    size: int, optional, default=None
+        size of the diagonal matrix (size x size).
+        Defaults to the number of spatial dimensions when unspecified
     """
     dim = size or len(func.dimensions)
     dim = dim-1 if func.is_TimeDependent else dim