Merge pull request #2768 from devitocodes/stagg-param

Stagg param

Author: mloubout

devito/finite_differences/differentiable.py

@@ -1133,7 +1133,9 @@ def _(expr, x0, **kwargs):
 def _(expr, x0, **kwargs):
     from devito.finite_differences.derivative import Derivative
     x0_expr = {d: v for d, v in x0.items() if v is not expr.indices_ref[d]}
-    if x0_expr:
+    if expr.is_parameter:
+        return expr
+    elif x0_expr:
         dims = tuple((d, 0) for d in x0_expr)
         fd_o = tuple([expr.interp_order]*len(dims))
         return Derivative(expr, *dims, fd_order=fd_o, x0=x0_expr)

devito/types/dense.py

@@ -301,7 +301,8 @@ def shape_global(self):
 
     @property
     def symbolic_shape(self):
-        return tuple(self._C_get_field(FULL, d).size for d in self.dimensions)
+        return DimensionTuple(*[self._C_get_field(FULL, d).size for d in self.dimensions],
+                              getters=self.dimensions)
 
     @property
     def size_global(self):
@@ -1018,7 +1019,7 @@ class Function(DiscreteFunction):
     is_autopaddable = True
 
     __rkwargs__ = (DiscreteFunction.__rkwargs__ +
-                   ('space_order', 'interp_order', 'dimensions'))
+                   ('space_order', 'interp_order', 'dimensions', 'is_parameter'))
 
     def _cache_meta(self):
         # Attach additional metadata to self's cache entry
@@ -1063,7 +1064,7 @@ def __init_finalize__(self, *args, **kwargs):
         # Used at operator evaluation to evaluate the Function at the
         # variable location (i.e. if the variable is staggered in x the
         # parameter has to be computed at x + hx/2)
-        self._is_parameter = kwargs.get('parameter', False)
+        self._is_parameter = kwargs.get('parameter', kwargs.get('is_parameter', False))
 
     def __fd_setup__(self):
         """

tests/test_derivatives.py

@@ -762,6 +762,24 @@ def test_nested_call(self):
         # Should be commutative
         assert simplify(deriv.evaluate - derivc.evaluate) == 0
 
+    def test_param_stagg_inner(self):
+        space_order = 2
+        nx, ny = 5, 5
+
+        grid = Grid((nx, ny))
+
+        x, y = grid.dimensions
+        yp = y + y.spacing / 2
+        xp = x + x.spacing / 2
+
+        f = TimeFunction(name="f", grid=grid, space_order=space_order, staggered=y)
+        p = Function(name="p", grid=grid, space_order=space_order, parameter=True)
+        g = TimeFunction(name="g", grid=grid, space_order=space_order, staggered=(x, y))
+
+        eqn = Eq(g, (p * f).dx)
+        eqne = eqn.evaluate.rhs
+        assert simplify(eqne - (p._subs(y, yp).evaluate * f).dx(x0=xp).evaluate) == 0
+
 
 class TestTwoStageEvaluation: