api: move symbolic config to its own options

Author: mloubout

devito/core/cpu.py

@@ -87,7 +87,6 @@ def _normalize_kwargs(cls, **kwargs):
 
         # Code generation options for derivatives
         o['expand'] = oo.pop('expand', cls.EXPAND)
-        o['interp-mode'] = oo.pop('interp-mode', cls.INTERP_MODE)
         o['deriv-collect'] = oo.pop('deriv-collect', cls.DERIV_COLLECT)
         o['deriv-schedule'] = oo.pop('deriv-schedule', cls.DERIV_SCHEDULE)
         o['deriv-unroll'] = oo.pop('deriv-unroll', False)
@@ -112,6 +111,7 @@ def _normalize_kwargs(cls, **kwargs):
             )
 
         kwargs['options'].update(o)
+        kwargs['sym_options'] = cls._normalize_sym_kwargs(**kwargs)
 
         return kwargs
 

devito/core/gpu.py

@@ -102,7 +102,6 @@ def _normalize_kwargs(cls, **kwargs):
 
         # Code generation options for derivatives
         o['expand'] = oo.pop('expand', cls.EXPAND)
-        o['interp-mode'] = oo.pop('interp-mode', cls.INTERP_MODE)
         o['deriv-collect'] = oo.pop('deriv-collect', cls.DERIV_COLLECT)
         o['deriv-schedule'] = oo.pop('deriv-schedule', cls.DERIV_SCHEDULE)
         o['deriv-unroll'] = oo.pop('deriv-unroll', False)
@@ -122,6 +121,7 @@ def _normalize_kwargs(cls, **kwargs):
             )
 
         kwargs['options'].update(o)
+        kwargs['sym_options'] = cls._normalize_sym_kwargs(**kwargs)
 
         return kwargs
 

devito/core/operator.py

@@ -125,26 +125,6 @@ class BasicOperator(Operator):
     finite-difference derivatives.
     """
 
-    INTERP_MODE = 'direct'
-    """
-    Interpolation mode used by `Mul._eval_at` when projecting a multi-factor
-    expression onto a target staggered location:
-
-    * `'direct'` (default): each factor is shifted to `func`'s location
-      independently (`Function._eval_at` per arg). Cheapest stencil; the
-      mode to pick unless you need an explicitly self-adjoint discretization.
-
-    * `'symmetric'`: when every factor lives at a staggering different from
-      `func`'s, the symmetric form `I * (a * I^T * b)` is built -- all
-      factors are gathered at the highest-priority "block" location via
-      `I^T`, multiplied there, and the product is interpolated to `func`
-      via `I`. Use this for operators whose continuous form decomposes as
-      `I * A * I^T` (e.g. the elastic stiffness `σ = C ε`).
-
-    See `examples/userapi/08_staggered_interp.ipynb` for the maths and a
-    worked elastic-stiffness example.
-    """
-
     DERIV_COLLECT = True
     """
     Factorize finite-difference derivatives exploiting the linearity of the FD
@@ -191,6 +171,30 @@ class BasicOperator(Operator):
     The target language constructor, to be specified by subclasses.
     """
 
+    # ------------------------------------------------------------------
+    # Symbolic-level option defaults (`sym_opt`).
+    # These steer mathematical choices made during expression lowering,
+    # *not* code generation or performance. They are kept separate from
+    # the `opt` options above to keep the two concerns distinct.
+    # ------------------------------------------------------------------
+
+    INTERP_MODE = 'direct'
+    """
+    Default for the `sym_opt={'interp-mode': ...}` option. Controls how
+    a product of fields living at different staggered locations is mapped
+    onto a target location:
+
+    * `'direct'` (default): each factor is interpolated to the target
+      independently. Cheapest stencil.
+    * `'symmetric'`: factors are first gathered at a common "block"
+      location, multiplied there, and the result is interpolated once to
+      the target. Preserves the `I A I^T` matrix structure, so the
+      discrete operator stays self-adjoint when the continuous one is
+      (e.g. the elastic stiffness `sigma = C eps`).
+
+    See `examples/userapi/08_staggered_interp.ipynb` for a worked example.
+    """
+
     @classmethod
     def _normalize_kwargs(cls, **kwargs):
         # Will be populated with dummy values; this method is actually overridden
@@ -208,12 +212,30 @@ def _normalize_kwargs(cls, **kwargs):
             )
 
         kwargs['options'].update(o)
+        kwargs['sym_options'] = cls._normalize_sym_kwargs(**kwargs)
 
         return kwargs
 
+    @classmethod
+    def _normalize_sym_kwargs(cls, **kwargs):
+        """
+        Fill in defaults and validate keys for the `sym_opt` dict passed to
+        the Operator. Returns the normalized `sym_options` dict.
+        """
+        so = dict(kwargs.get('sym_options', {}))
+        out = {'interp-mode': so.pop('interp-mode', cls.INTERP_MODE)}
+
+        if so:
+            raise InvalidOperator(
+                f'Unrecognized symbolic options: [{", ".join(list(so))}]'
+            )
+
+        return out
+
     @classmethod
     def _check_kwargs(cls, **kwargs):
         oo = kwargs['options']
+        so = kwargs['sym_options']
 
         if oo['mpi'] and oo['mpi'] not in cls.MPI_MODES:
             raise InvalidOperator(f"Unsupported MPI mode `{oo['mpi']}`")
@@ -229,6 +251,9 @@ def _check_kwargs(cls, **kwargs):
         if oo['errctl'] not in (None, False, 'basic', 'max'):
             raise InvalidOperator("Illegal `errctl` value")
 
+        if so['interp-mode'] not in ('direct', 'symmetric'):
+            raise InvalidOperator("Illegal `interp-mode` value")
+
     def _autotune(self, args, setup):
         if setup in [False, 'off']:
             return args

devito/finite_differences/differentiable.py

@@ -1,5 +1,4 @@
 from collections import ChainMap
-from contextlib import suppress
 from functools import cached_property, singledispatch
 from itertools import product
 
@@ -16,6 +15,7 @@
     # Moved in 1.13
     from sympy.core.basic import ordering_of_classes
 
+from devito.finite_differences.interpolation import interp_at, post_x0_indices
 from devito.finite_differences.tools import coeff_priority, make_shift_x0
 from devito.logger import warning
 from devito.tools import (
@@ -699,32 +699,31 @@ def _eval_at(self, func, interp_mode='direct', **kwargs):
         if interp_mode != 'symmetric':
             return super()._eval_at(func, **kwargs)
 
-        diff_args = [a for a in self.args if isinstance(a, Differentiable)]
-        other_args = [a for a in self.args if not isinstance(a, Differentiable)]
+        diff, other = split(self.args, lambda a: isinstance(a, Differentiable))
 
         # Symmetric form requires every Differentiable factor to differ from
         # func; otherwise direct evaluation is cleaner and equivalent.
-        if len(diff_args) < 2 or \
-           any(a.staggered == func.staggered for a in diff_args):
+        if len(diff) < 2 or \
+           any(a.staggered == func.staggered for a in diff):
             return super()._eval_at(func, **kwargs)
 
         block_indices = highest_priority(self).indices_ref
 
         # Bring each factor to block's location (I^T where needed)
-        new_factors = list(other_args)
-        for a in diff_args:
+        new_factors = list(other)
+        for a in diff:
             if isinstance(a, sympy.Derivative):
-                source = _post_x0_indices(a, func)
+                source = post_x0_indices(a, func)
                 a = a._rebuild(x0={dim: func.indices_ref[dim] for dim in a.dims
                                    if dim in func.indices_ref.getters})
             else:
                 source = a.indices_ref
-            new_factors.append(_interp_at(a, source, block_indices,
-                                          self.interp_order))
+            new_factors.append(interp_at(a, source, block_indices,
+                                         self.interp_order))
 
         # Final I from block's location to func
-        return _interp_at(self.func(*new_factors), block_indices,
-                          func.indices_ref, self.interp_order)
+        return interp_at(self.func(*new_factors), block_indices,
+                         func.indices_ref, self.interp_order)
 
 
 class Pow(DifferentiableOp, sympy.Pow):
@@ -1251,63 +1250,6 @@ def _diff2sympy(obj):
 evalf_table[Pow] = evalf_table[sympy.Pow]
 
 
-def _interp_mapper(source, target, dims):
-    """
-    Build a `{dim: target_index}` mapper for dimensions in `dims` where
-    `source[dim]` differs from `target[dim]`.
-
-    `source` and `target` are dict-like `{dim: index_expr}` (e.g. a plain
-    dict or a `DimensionTuple`). Dimensions missing from either side are
-    skipped silently.
-    """
-    mapper = {}
-    for d in dims:
-        try:
-            s = source[d]
-            t = target[d]
-        except (KeyError, IndexError):
-            continue
-        if s is not t:
-            mapper[d] = t
-    return mapper
-
-
-def _interp_at(expr, source, target, interp_order):
-    """
-    Build a symbolic 0-order FD interpolation operator on `expr` that maps
-    values from `source` indices to `target` indices, only on the
-    dimensions where the two locations differ.
-    """
-    if not isinstance(expr, Differentiable):
-        return expr
-
-    mapper = _interp_mapper(source, target, expr.dimensions)
-    if not mapper:
-        return expr
-
-    return expr.diff(*mapper.keys(),
-                     deriv_order=(0,) * len(mapper),
-                     fd_order=(interp_order,) * len(mapper),
-                     x0=mapper)
-
-
-def _post_x0_indices(deriv, func):
-    """
-    Conceptual indices of `deriv` after setting `x0` on its own derivative
-    dimensions to `func`'s indices. Derivative dims take `func`'s indices;
-    other dims keep the underlying expression's natural location (so that
-    `interp_for_fd` does not introduce a spurious second shift).
-    """
-    ref = {}
-    for dim in deriv.dimensions:
-        if dim in deriv.dims and dim in func.indices_ref.getters:
-            ref[dim] = func.indices_ref[dim]
-        else:
-            with suppress(KeyError):
-                ref[dim] = deriv.indices_ref[dim]
-    return ref
-
-
 # Interpolation for finite differences
 @singledispatch
 def interp_for_fd(expr, x0, **kwargs):

devito/finite_differences/interpolation.py

@@ -0,0 +1,62 @@
+from contextlib import suppress
+
+__all__ = ['interp_at', 'interp_mapper', 'post_x0_indices']
+
+
+def interp_mapper(source, target, dims):
+    """
+    Build a `{dim: target_index}` mapper for dimensions in `dims` where
+    `source[dim]` differs from `target[dim]`.
+
+    `source` and `target` are dict-like `{dim: index_expr}` (e.g. a plain
+    dict or a `DimensionTuple`). Dimensions missing from either side are
+    skipped silently.
+    """
+    mapper = {}
+    for d in dims:
+        try:
+            s = source[d]
+            t = target[d]
+        except (KeyError, IndexError):
+            continue
+        if s is not t:
+            mapper[d] = t
+    return mapper
+
+
+def interp_at(expr, source, target, interp_order):
+    """
+    Build a symbolic 0-order FD interpolation operator on `expr` that maps
+    values from `source` indices to `target` indices, only on the
+    dimensions where the two locations differ.
+    """
+    from devito.finite_differences.differentiable import Differentiable
+
+    if not isinstance(expr, Differentiable):
+        return expr
+
+    mapper = interp_mapper(source, target, expr.dimensions)
+    if not mapper:
+        return expr
+
+    return expr.diff(*mapper.keys(),
+                     deriv_order=(0,) * len(mapper),
+                     fd_order=(interp_order,) * len(mapper),
+                     x0=mapper)
+
+
+def post_x0_indices(deriv, func):
+    """
+    Conceptual indices of `deriv` after setting `x0` on its own derivative
+    dimensions to `func`'s indices. Derivative dims take `func`'s indices;
+    other dims keep the underlying expression's natural location (so that
+    `interp_for_fd` does not introduce a spurious second shift).
+    """
+    ref = {}
+    for dim in deriv.dimensions:
+        if dim in deriv.dims and dim in func.indices_ref.getters:
+            ref[dim] = func.indices_ref[dim]
+        else:
+            with suppress(KeyError):
+                ref[dim] = deriv.indices_ref[dim]
+    return ref

devito/operator/operator.py

@@ -67,6 +67,16 @@ class Operator(Callable):
             Symbolic substitutions to be applied to ``expressions``.
         * opt : str
             The performance optimization level. Defaults to ``configuration['opt']``.
+        * sym_opt : dict
+            Symbolic-level options controlling mathematical choices made during
+            expression lowering (e.g. how staggered multi-factor products are
+            interpolated). Distinct from ``opt``, which controls code generation
+            and performance. Accepted keys:
+
+            - ``'interp-mode'`` (``'direct'`` | ``'symmetric'``): selects the
+              interpolation strategy used by ``Mul._eval_at`` when projecting a
+              multi-factor expression onto a target staggered location. See the
+              tutorial at ``examples/userapi/08_staggered_interp.ipynb``.
         * language : str
             The target language for shared-memory parallelism. Defaults to
             ``configuration['language']``.
@@ -234,6 +244,7 @@ def _build(cls, expressions, **kwargs):
         # Potentially required for lazily allocated Functions
         op._mode = kwargs['mode']
         op._options = kwargs['options']
+        op._sym_options = kwargs['sym_options']
         op._allocator = kwargs['allocator']
         op._platform = kwargs['platform']
 
@@ -341,7 +352,7 @@ def _lower_exprs(cls, expressions, **kwargs):
             * Shift indices for domain alignment.
         """
         expand = kwargs['options'].get('expand', True)
-        interp_mode = kwargs['options'].get('interp-mode', 'direct')
+        interp_mode = kwargs.get('sym_options', {}).get('interp-mode', 'direct')
 
         # Specialization is performed on unevaluated expressions
         expressions = cls._specialize_dsl(expressions, **kwargs)
@@ -1643,6 +1654,12 @@ def parse_kwargs(**kwargs):
         mode = 'noop'
     kwargs['mode'] = mode
 
+    # `sym_opt` -- symbolic-level options (mathematical choices, not codegen)
+    sym_opt = kwargs.pop('sym_opt', None) or {}
+    if not isinstance(sym_opt, (dict, frozendict)):
+        raise InvalidOperator(f"Illegal `sym_opt={str(sym_opt)}`")
+    kwargs['sym_options'] = dict(sym_opt)
+
     # `platform`
     platform = kwargs.get('platform')
     if platform is not None:

devito/types/dense.py

@@ -15,7 +15,7 @@
 from devito.deprecations import deprecations
 from devito.exceptions import InvalidArgument
 from devito.finite_differences import Differentiable, generate_fd_shortcuts
-from devito.finite_differences.differentiable import _interp_mapper
+from devito.finite_differences.interpolation import interp_mapper
 from devito.finite_differences.tools import fd_weights_registry
 from devito.logger import debug, warning
 from devito.mpi import MPI
@@ -1128,7 +1128,7 @@ def _eval_at(self, func, **kwargs):
             return self
 
         # Dims where self and func indices differ -> {dim: func_idx}
-        diff = _interp_mapper(self.indices_ref, func.indices_ref, self.dimensions)
+        diff = interp_mapper(self.indices_ref, func.indices_ref, self.dimensions)
 
         # Translate into a subs mapper {self_idx: func_idx} aligned on self's dims
         subs_map = {self.indices_ref[d]: t._subs(func.dimensions[d], d)