From d5938784ab89f659b2422c43164568c16997ad62 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 21 Oct 2026 11:52:27 -0400 Subject: [PATCH] compiler: Interpolate non-integer SubGrid reads as a plain stencil A SubGrid Function used in an equation is now read at the left-hand side's location, the convention staggered Functions already follow -- no wrapper and no flag: sub = SubGrid(grid, factor=1.5) vp = Function(name='vp', grid=sub) Eq(u, vp*u) `SubGrid` accepts non-integer and sub-unit factors, dispatching on the factor: integer factors keep the strided ConditionalDimension lowering, anything else gets its own SpaceDimensions and interpolates. `is_interpolated` distinguishes the two, and a factor of exactly 1 passes the parent's Dimension through. `Function._eval_at` emits a `SubGridInterp` marker that FD lowering leaves inert, and the new Cluster pass `subgrid_interp` expands it once the iteration space is known -- the same deferral `IndexDerivative` uses. Reading a SubGrid at the parent's point `x` evaluates at `x*q/p` for `factor = p/q`, so the base index is `(x*q)//p` and the fractional part repeats with period `p`. The weights are therefore constants within a phase and live in a compile-time table indexed by `Mod(x, p)`: no phase loop, no branch, and unit-stride loops over the parent's Dimensions alone. --- devito/finite_differences/interpolation.py | 100 ++- devito/operator/operator.py | 10 +- devito/passes/clusters/__init__.py | 1 + devito/passes/clusters/subgrid.py | 193 ++++ devito/types/dense.py | 13 +- devito/types/grid.py | 111 ++- examples/userapi/10_subgrids.ipynb | 967 +++++++++++++++++---- tests/test_subgrid.py | 288 ++++++ tests/test_subsampling.py | 24 +- 9 files changed, 1463 insertions(+), 244 deletions(-) create mode 100644 devito/passes/clusters/subgrid.py create mode 100644 tests/test_subgrid.py diff --git a/devito/finite_differences/interpolation.py b/devito/finite_differences/interpolation.py index 7fb397b080..74e8011686 100644 --- a/devito/finite_differences/interpolation.py +++ b/devito/finite_differences/interpolation.py @@ -1,6 +1,104 @@ from contextlib import suppress -__all__ = ['interp_at', 'interp_mapper', 'post_x0_indices'] +import sympy + +from devito.types.lazy import Evaluable + +__all__ = ['SubGridInterp', 'interp_at', 'interp_mapper', 'is_subgrid_read', + 'post_x0_indices'] + + +class SubGridInterp(sympy.Expr, Evaluable): + + """ + A SubGrid Function read at the location of a parent-Grid Function, + unevaluated. + + A SubGrid point `i` sits at parent coordinate `i*factor`, so reading it at + the parent point `x` means evaluating at `x/factor` -- a fractional offset + that varies with the iteration index, which the finite-difference + machinery cannot turn into a compile-time stencil window. + + What makes it a stencil after all is that the offset is *periodic*. With + `factor = p/q` the coordinate is `x*q/p`, so + + base = (x*q) // p frac = ((x*q) % p) / p + + and `frac` depends on `x` only through `x % p`: there are `p` phases, and + within a phase the tap weights are constants. Deriving that needs the + iteration space, so it is done by a Cluster pass + (`passes/clusters/subgrid.py`); this object just records *what* to relocate + and survives FD lowering untouched, in the manner of `IndexDerivative`. + + Parameters + ---------- + func : Function + The SubGrid Function being read. + target : Function + The Function whose location to read it at. + """ + + __rargs__ = ('func', 'target') + + is_commutative = True + + def __new__(cls, func, target=None, **kwargs): + # `target` rides along as a sympy arg so that generic reconstruction + # (`expr.func(*expr.args)`, as `uxreplace` does) round-trips + obj = sympy.Expr.__new__(cls, func, target) + obj._func = func + obj._target = target + + return obj + + @property + def func_(self): + """The SubGrid Function being read.""" + return self._func + + @property + def target(self): + """The Function whose location it is read at.""" + return self._target + + @property + def subgrid(self): + return self._func.function.grid + + @property + def dtype(self): + return self._func.function.dtype + + def __repr__(self): + return f"SubGridInterp({self._func}, at={self._target})" + + __str__ = __repr__ + + def _evaluate(self, **kwargs): + # Deliberately inert: the relocation restructures the iteration space, + # so it is lowered by a Cluster pass + return self + + +def is_subgrid_read(func, target): + """ + Whether reading `func` at `target`'s location is a SubGrid relocation. + + True when `func` lives on a SubGrid and `target` does not live on that same + SubGrid. The exception matters: a Function read at a shift of its own + Dimensions, as a recurrence does, is not a relocation. + + Only the Dimensions `target` actually iterates can be relocated, so a + SubGrid sharing none of them is not a relocation either. + """ + grid = getattr(func, 'grid', None) + if grid is None or not getattr(grid, 'is_SubGrid', False): + return False + + if getattr(target, 'grid', None) is grid: + return False + + return any(d in target.dimensions for d in grid.parent.dimensions) def interp_mapper(source, target, dims): diff --git a/devito/operator/operator.py b/devito/operator/operator.py index eb6a17326e..19880a306d 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -34,7 +34,7 @@ from devito.passes import ( Graph, error_mapper, finalize_args, generate_implicit, generate_macros, is_on_device, lower_dtypes, lower_index_derivatives, lower_local_sums, minimize_symbols, - optimize_pows, subsampling, unevaluate + optimize_pows, subgrid_interp, subsampling, unevaluate ) from devito.symbolics import estimate_cost, subs_op_args from devito.tools import ( @@ -423,9 +423,11 @@ def _lower_clusters(cls, expressions, profiler=None, **kwargs): # Generate implicit Clusters from higher level abstractions clusters = generate_implicit(clusters) - # Strengthen factor-based ConditionalDimension guards into strided - # iteration. Before `lower_index_derivatives` so that the Dimension - # swap is seen by everything downstream + # Lower SubGrid reads into stencils, and strengthen factor-based + # ConditionalDimension guards into strided iteration. Both rewrite the + # iteration space, so they run before `lower_index_derivatives` and + # everything else downstream sees the result + clusters = subgrid_interp(clusters, **kwargs) clusters = subsampling(clusters) # Lower all remaining high order symbolic objects diff --git a/devito/passes/clusters/__init__.py b/devito/passes/clusters/__init__.py index 3a6cb0d060..7c6311fcd0 100644 --- a/devito/passes/clusters/__init__.py +++ b/devito/passes/clusters/__init__.py @@ -6,6 +6,7 @@ from .fusion import * # noqa from .blocking import * # noqa from .subsampling import * # noqa +from .subgrid import * # noqa from .asynchrony import * # noqa from .implicit import * # noqa from .misc import * # noqa diff --git a/devito/passes/clusters/subgrid.py b/devito/passes/clusters/subgrid.py new file mode 100644 index 0000000000..236eb36d81 --- /dev/null +++ b/devito/passes/clusters/subgrid.py @@ -0,0 +1,193 @@ +from fractions import Fraction + +import sympy + +from devito.finite_differences.interpolation import SubGridInterp +from devito.ir.clusters import Queue +from devito.symbolics import search, uxreplace +from devito.tools import timed_pass +from devito.types import Array, CustomDimension + +__all__ = ['subgrid_interp'] + + +@timed_pass() +def subgrid_interp(clusters, sregistry=None, **kwargs): + """ + Lower a Cluster reading an interpolated SubGrid Function into a plain + stencil. + + A SubGrid point `i` sits at parent coordinate `i*factor`, so reading it at + the parent point `x` means evaluating at `x/factor`. With `factor = p/q` + that coordinate is `x*q/p`, hence + + base = (x*q) // p frac = ((x*q) % p) / p + + Both are plain integer arithmetic, and `frac` depends on `x` only through + `x % p` -- so there are `p` phases, and within a phase the tap weights are + *constants*. Only the weights vary between phases, so they go in a small + compile-time table indexed by `Mod(x, p)`. + + No phase loop is therefore needed, and nothing is added to the iteration + space: the SubGrid Dimensions are simply dropped and the inner loops stay + unit-stride and branch-free. + """ + return SubGridInterpolation(sregistry).process(clusters) + + +class SubGridInterpolation(Queue): + + def __init__(self, sregistry=None, **kwargs): + super().__init__(**kwargs) + self.sregistry = sregistry + + def callback(self, clusters, prefix): + if prefix: + return clusters + + processed = [] + for c in clusters: + processed.append(_lower(c, self.sregistry)) + + return processed + + +def _lower(cluster, sregistry): + """Expand the SubGridInterp markers in one Cluster, or return it as is.""" + markers = set() + for e in cluster.exprs: + markers.update(search(e, SubGridInterp)) + + if not markers: + return cluster + + ispace = cluster.ispace + weights = {} # (factor, dtype) -> (Array, ntaps, nphases) + subs = {} + + # Blocking may already have run, in which case a parent Dimension appears + # in the IterationSpace only through its block Dimensions (`y -> y0_blk0`). + # Promote those back so the membership test below sees the roots + itdims = set(ispace.promote(lambda d: d.is_Block).itdims) + + for m in markers: + grid = m.subgrid + func = m.func_.function + + terms = [] + for ds, d, fac in zip(grid.dimensions, grid.parent.dimensions, + grid.factor, strict=True): + if d not in itdims: + continue + + frac = Fraction(float(fac)).limit_denominator(10**6) + key = (float(fac), func.dtype) + if key not in weights: + name = sregistry.make_name(prefix='w') if sregistry else f'w{d.name}' + weights[key] = _phase_weights(fac, func.dtype, name) + w, ntaps, nphases = weights[key] + + # `nphases == 1` means every parent point lands exactly on a + # SubGrid point: a single tap of unit weight, no interpolation + ph = sympy.Mod(d, nphases) if nphases > 1 else None + terms.append((ds, d, ph, w, ntaps, frac)) + + if not terms: + continue + + # The relocated read: a weighted sum over the taps of every relocated + # Dimension, with the weights taken from each Dimension's phase row + subs[m] = _stencil(func, terms) + + exprs = [uxreplace(e, subs) for e in cluster.exprs] + + # The SubGrid Dimensions are no longer iterated: the relocated reads are + # indexed by the parent's Dimensions alone. An axis of factor 1 is the + # exception -- there the SubGrid shares the parent's own Dimension, which + # stays in the IterationSpace + gone = {ds for m in markers + for ds, d in zip(m.subgrid.dimensions, m.subgrid.parent.dimensions, + strict=True) + if ds is not d} + ispace = ispace.project(lambda d: d not in gone) + + return cluster.rebuild(exprs=exprs, ispace=ispace) + + +def _stencil(func, terms): + """ + `func` read as a weighted sum over the taps of every relocated Dimension, + with the weights taken from each Dimension's phase row. + """ + from itertools import product + + retval = [] + ranges = [range(t[4]) if t[2] is not None else range(1) for t in terms] + for offsets in product(*ranges): + weight = sympy.S.One + indices = {} + for (ds, d, ph, w, ntaps, frac), o in zip(terms, offsets, + strict=True): + if ph is not None: + weight = weight*w[ntaps*ph + o] + indices[ds] = _base(d, frac) + o + + # `func` is indexed directly here, which bypasses the halo shift that + # index lowering applies to ordinary reads -- so add it explicitly + idxs = [indices.get(d, d) + h for d, h in + zip(func.dimensions, func._size_nodomain.left, strict=True)] + retval.append(weight*func.indexed[tuple(idxs)]) + + return sympy.Add(*retval) + + +def _base(d, frac): + """ + The SubGrid base index for parent index `d`: `(d*q)//p` for a + `factor = p/q`, in integer arithmetic. + + `d` is the raw iteration Dimension: index lowering adds each Function's + halo afterwards, so neither halo takes part in the division. + """ + from devito.symbolics import IntDiv + return IntDiv(frac.denominator*d, frac.numerator) + + +def _phase_weights(factor, dtype, name): + """ + The per-phase tap weights for *reading* a SubGrid of the given `factor` at + its parent's points, as a compile-time table. + + A SubGrid point sits at parent coordinate `i*factor`, so conversely the + parent point `x` sits at SubGrid coordinate `x/factor`. With + `factor = p/q` that is `x*q/p`, hence + + base = (x*q)//p frac = ((x*q) % p)/p + + `base` is the integer part, and `frac` depends on `x` only through + `x % p`. So there are `p` phases, and within a phase the two linear + weights `(1 - frac, frac)` are constants -- an ordinary constant-offset + stencil, read at `base` and `base + 1`. + + Returns `(weights, ntaps, nphases)`: the `Array` holding `nphases*ntaps` + values, the taps per phase, and the number of phases. + """ + frac = Fraction(float(factor)).limit_denominator(10**6) + p, q = frac.numerator, frac.denominator + + # Reading the SubGrid at the parent's points inverts the factor, so the + # phase period is `p` + nphases, ntaps = p, 2 + + values = [] + for ph in range(nphases): + f = float(Fraction((ph*q) % p, p)) + values.extend(sympy.Float(v) for v in (1. - f, f)) + + d = CustomDimension(name=f'{name}d', symbolic_min=0, + symbolic_max=nphases*ntaps - 1, + symbolic_size=nphases*ntaps) + weights = Array(name=name, dimensions=(d,), shape=(nphases*ntaps,), + dtype=dtype, scope='stack', initvalue=values) + + return weights, ntaps, nphases diff --git a/devito/types/dense.py b/devito/types/dense.py index 30f543a47c..ad7b5d838a 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -15,7 +15,9 @@ from devito.deprecations import deprecations from devito.exceptions import InvalidArgument from devito.finite_differences import Differentiable, generate_fd_shortcuts -from devito.finite_differences.interpolation import interp_mapper +from devito.finite_differences.interpolation import ( + SubGridInterp, interp_mapper, is_subgrid_read +) from devito.finite_differences.tools import fd_weights_registry from devito.logger import debug, warning from devito.mpi import MPI @@ -1130,6 +1132,15 @@ def _fd_priority(self): return 1.2 if self.staggered.on_node else 1.1 def _eval_at(self, func, **kwargs): + # A Function on a SubGrid lives on different Dimensions altogether, so + # relocating it onto `func` means re-expressing its indices in the + # parent's terms. Same convention as staggering: the RHS is evaluated + # at the LHS's location + if is_subgrid_read(self, func): + # The relocation depends on the iteration index, so it cannot be + # turned into a stencil here; defer it to the Cluster pass + return SubGridInterp(self, func) + if self.staggered == func.staggered or self.interp_order == 0: return self diff --git a/devito/types/grid.py b/devito/types/grid.py index 4568b852bf..cb4e7f680d 100644 --- a/devito/types/grid.py +++ b/devito/types/grid.py @@ -12,7 +12,7 @@ from devito.finite_differences.elementary import Max, Min from devito.logger import warning from devito.mpi import MPI, Distributor, SubDistributor -from devito.tools import ReducerMap, as_tuple, frozendict, is_integer +from devito.tools import ReducerMap, as_tuple, frozendict, is_integer, is_number from devito.types.args import ArgProvider from devito.types.basic import Scalar from devito.types.dense import Function @@ -36,6 +36,7 @@ class CartesianDiscretization: """ is_Grid = False + is_SubGrid = False is_SubDomain = False is_MultiSubDomain = False @@ -432,23 +433,28 @@ def __setstate__(self, state): class SubGrid(Grid): """ - A coarser Grid whose points are a subset of a parent Grid's points. + A Grid whose points are derived from a parent Grid's by a factor. - Along each Dimension subsampled by `factor`, the SubGrid retains every - `factor`-th point of the parent, so point `i` of the SubGrid sits at parent - index `i*factor`. The physical spacing is `factor` times the parent's, and - the origin is shared. + Along each Dimension scaled by `factor`, point `i` of the SubGrid sits at + parent index `i*factor`. The physical spacing is `factor` times the + parent's, and the origin is shared. + + An **integer** factor puts SubGrid points exactly on parent points, so the + SubGrid is a subsample of the parent and the compiler reaches it by + striding. A **non-integer** factor puts them in between, so reading such a + Function on the parent grid interpolates; `is_interpolated` distinguishes + the two. A factor below 1 refines rather than coarsens. Parameters ---------- parent : Grid - The Grid being subsampled. - factor : int or tuple of int - Subsampling factor, either one value per space Dimension of `parent`, - or a single value applied to all of them. A factor of 1 leaves the - corresponding Dimension untouched. + The Grid being rescaled. + factor : int or float or tuple + Scaling factor, either one value per space Dimension of `parent`, or a + single value applied to all of them. Must be positive. A factor of 1 + leaves the corresponding Dimension untouched. name : str, optional - Prefix for the names of the subsampled Dimensions; each is named + Prefix for the names of the derived Dimensions; each is named `f'{parent_dimension}{name}'`. Defaults to `'s'`, giving `xs`, `ys`. Examples @@ -483,11 +489,18 @@ class SubGrid(Grid): Notes ----- - Because the subsampled points coincide with parent points, no interpolation - is involved: an `Eq` writing a SubGrid Function from a parent Function is a - strided copy, and the compiler emits a strided loop for it. + With an integer factor the points coincide with parent points, so no + interpolation is involved: an `Eq` writing a SubGrid Function from a parent + Function is a strided copy, and the compiler emits a strided loop for it. + + With a non-integer factor, a SubGrid Function appearing in an equation is + read at the *left-hand side's* location, the same convention a staggered + Function follows. It lowers to a plain stencil whose weights come from a + compile-time table indexed by the phase, so the loops stay unit-stride. """ + is_SubGrid = True + def __init__(self, parent, factor, name='s', **kwargs): if not parent.is_Grid: raise ValueError("`parent` must be a Grid, got " @@ -508,17 +521,31 @@ def __init__(self, parent, factor, name='s', **kwargs): f"`SubGrid(parent.parent, factor={combined})`" ) - dimensions = tuple( - d if f == 1 else ConditionalDimension(f'{d.name}{name}', parent=d, - factor=f) - for d, f in zip(parent.dimensions, factors, strict=True) - ) - - # A subsampled axis keeps the parent's endpoints, so the number of - # retained points is `ceil(n/factor)` and the extent is unchanged up to - # the last retained point - shape = tuple((n + f - 1)//f - for n, f in zip(parent.shape, factors, strict=True)) + if all(is_integer(f) for f in factors): + # Points coincide with the parent's: derive ConditionalDimensions, + # so a plain `Eq` lowers to a strided copy + dimensions = tuple( + d if f == 1 else ConditionalDimension(f'{d.name}{name}', + parent=d, factor=f) + for d, f in zip(parent.dimensions, factors, strict=True) + ) + # A subsampled axis keeps the parent's endpoints, so exactly + # `ceil(n/factor)` points are retained + shape = tuple((n + f - 1)//f + for n, f in zip(parent.shape, factors, strict=True)) + else: + # Points fall between the parent's, so the axis is a Dimension in + # its own right and reading it on the parent interpolates + dimensions = tuple( + d if f == 1 else SpaceDimension( + f'{d.name}{name}', + spacing=Spacing(name=f'h_{d.name}{name}', is_const=True) + ) + for d, f in zip(parent.dimensions, factors, strict=True) + ) + # The last point must not overshoot the parent's domain + shape = tuple(int(math.floor((n - 1)/f)) + 1 + for n, f in zip(parent.shape, factors, strict=True)) extent = tuple(s*(n - 1) for s, n in zip(self.__spacing_setup__(parent, factors), shape, strict=True)) @@ -548,22 +575,16 @@ def __factors_setup__(cls, parent, factor): f"got {len(factors)}") for f in factors: - if not is_integer(f): - # Non-integer factors place SubGrid points between parent - # points, which requires interpolation rather than striding - raise ValueError( - f"Non-integer subsampling factor `{f}`. A SubGrid's points " - "must coincide with its parent's, so the factor must be a " - f"whole number of parent cells -- use " - f"`{max(1, math.floor(f))}` or `{max(1, math.ceil(f))}`. " - "Subsampling onto points that fall " - "between parent points requires interpolation, which " - "SubGrid does not perform" - ) - if f < 1: - raise ValueError(f"Subsampling factor must be >= 1, got `{f}`") + if not is_number(f): + raise ValueError(f"Subsampling factor must be a number, got `{f}`") + if f <= 0: + raise ValueError(f"Subsampling factor must be > 0, got `{f}`") - return tuple(int(f) for f in factors) + # Integer factors place SubGrid points exactly on parent points, which + # the compiler can reach by striding. Anything else falls between them + # and is read by interpolating + return tuple(int(f) if is_integer(f) or float(f).is_integer() + else float(f) for f in factors) @classmethod def __spacing_setup__(cls, parent, factors): @@ -588,6 +609,14 @@ def root(self): # under MPI return self._parent.root + @property + def is_interpolated(self): + """ + True if the SubGrid's points fall between its parent's, so that reading + such a Function on the parent grid interpolates rather than striding. + """ + return not all(is_integer(f) for f in self._factors) + def __repr__(self): return f"SubGrid[shape={self.shape}, factor={tuple(self._factors)}]" diff --git a/examples/userapi/10_subgrids.ipynb b/examples/userapi/10_subgrids.ipynb index 5a2b9b92df..3bdea57b9d 100644 --- a/examples/userapi/10_subgrids.ipynb +++ b/examples/userapi/10_subgrids.ipynb @@ -37,10 +37,10 @@ "execution_count": 1, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:09.864846Z", - "iopub.status.busy": "2026-10-21T00:00:09.864277Z", - "iopub.status.idle": "2026-10-21T00:00:10.707068Z", - "shell.execute_reply": "2026-10-21T00:00:10.705504Z" + "iopub.execute_input": "2026-10-21T15:11:18.753387Z", + "iopub.status.busy": "2026-10-21T15:11:18.752790Z", + "iopub.status.idle": "2026-10-21T15:11:19.504262Z", + "shell.execute_reply": "2026-10-21T15:11:19.503361Z" } }, "outputs": [ @@ -80,10 +80,10 @@ "execution_count": 2, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.761033Z", - "iopub.status.busy": "2026-10-21T00:00:10.760169Z", - "iopub.status.idle": "2026-10-21T00:00:10.769237Z", - "shell.execute_reply": "2026-10-21T00:00:10.767660Z" + "iopub.execute_input": "2026-10-21T15:11:19.545444Z", + "iopub.status.busy": "2026-10-21T15:11:19.545043Z", + "iopub.status.idle": "2026-10-21T15:11:19.550908Z", + "shell.execute_reply": "2026-10-21T15:11:19.549700Z" } }, "outputs": [ @@ -112,10 +112,10 @@ "execution_count": 3, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.775239Z", - "iopub.status.busy": "2026-10-21T00:00:10.774599Z", - "iopub.status.idle": "2026-10-21T00:00:10.787096Z", - "shell.execute_reply": "2026-10-21T00:00:10.785379Z" + "iopub.execute_input": "2026-10-21T15:11:19.556470Z", + "iopub.status.busy": "2026-10-21T15:11:19.555880Z", + "iopub.status.idle": "2026-10-21T15:11:19.568657Z", + "shell.execute_reply": "2026-10-21T15:11:19.567009Z" } }, "outputs": [ @@ -150,10 +150,10 @@ "execution_count": 4, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.792743Z", - "iopub.status.busy": "2026-10-21T00:00:10.792135Z", - "iopub.status.idle": "2026-10-21T00:00:10.802173Z", - "shell.execute_reply": "2026-10-21T00:00:10.800578Z" + "iopub.execute_input": "2026-10-21T15:11:19.574075Z", + "iopub.status.busy": "2026-10-21T15:11:19.573557Z", + "iopub.status.idle": "2026-10-21T15:11:19.582823Z", + "shell.execute_reply": "2026-10-21T15:11:19.581308Z" } }, "outputs": [ @@ -185,10 +185,10 @@ "execution_count": 5, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.807643Z", - "iopub.status.busy": "2026-10-21T00:00:10.807100Z", - "iopub.status.idle": "2026-10-21T00:00:10.817272Z", - "shell.execute_reply": "2026-10-21T00:00:10.815647Z" + "iopub.execute_input": "2026-10-21T15:11:19.588160Z", + "iopub.status.busy": "2026-10-21T15:11:19.587643Z", + "iopub.status.idle": "2026-10-21T15:11:19.597745Z", + "shell.execute_reply": "2026-10-21T15:11:19.596155Z" } }, "outputs": [ @@ -220,10 +220,10 @@ "execution_count": 6, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.822983Z", - "iopub.status.busy": "2026-10-21T00:00:10.822389Z", - "iopub.status.idle": "2026-10-21T00:00:10.836328Z", - "shell.execute_reply": "2026-10-21T00:00:10.834740Z" + "iopub.execute_input": "2026-10-21T15:11:19.603424Z", + "iopub.status.busy": "2026-10-21T15:11:19.602804Z", + "iopub.status.idle": "2026-10-21T15:11:19.616911Z", + "shell.execute_reply": "2026-10-21T15:11:19.615440Z" } }, "outputs": [ @@ -255,10 +255,10 @@ "execution_count": 7, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.841993Z", - "iopub.status.busy": "2026-10-21T00:00:10.841397Z", - "iopub.status.idle": "2026-10-21T00:00:10.852733Z", - "shell.execute_reply": "2026-10-21T00:00:10.851021Z" + "iopub.execute_input": "2026-10-21T15:11:19.622389Z", + "iopub.status.busy": "2026-10-21T15:11:19.621781Z", + "iopub.status.idle": "2026-10-21T15:11:19.632672Z", + "shell.execute_reply": "2026-10-21T15:11:19.631230Z" } }, "outputs": [ @@ -289,10 +289,10 @@ "execution_count": 8, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.858751Z", - "iopub.status.busy": "2026-10-21T00:00:10.858250Z", - "iopub.status.idle": "2026-10-21T00:00:10.867649Z", - "shell.execute_reply": "2026-10-21T00:00:10.866250Z" + "iopub.execute_input": "2026-10-21T15:11:19.637729Z", + "iopub.status.busy": "2026-10-21T15:11:19.637246Z", + "iopub.status.idle": "2026-10-21T15:11:19.646092Z", + "shell.execute_reply": "2026-10-21T15:11:19.644608Z" } }, "outputs": [ @@ -317,9 +317,10 @@ "source": [ "### What is rejected\n", "\n", - "Only factors that keep subgrid points on parent points are allowed. A non-integer\n", - "factor would place them in between, which needs interpolation rather than\n", - "subsampling." + "An integer factor keeps subgrid points *on* parent points, which is plain\n", + "subsampling. A non-integer factor places them in between; that is also supported,\n", + "through interpolation, and is covered further down. What is rejected is a factor\n", + "that is not positive, or the wrong number of them." ] }, { @@ -327,10 +328,10 @@ "execution_count": 9, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.872862Z", - "iopub.status.busy": "2026-10-21T00:00:10.872263Z", - "iopub.status.idle": "2026-10-21T00:00:10.882012Z", - "shell.execute_reply": "2026-10-21T00:00:10.880300Z" + "iopub.execute_input": "2026-10-21T15:11:19.651225Z", + "iopub.status.busy": "2026-10-21T15:11:19.650643Z", + "iopub.status.idle": "2026-10-21T15:11:19.659922Z", + "shell.execute_reply": "2026-10-21T15:11:19.658309Z" } }, "outputs": [ @@ -338,14 +339,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "factor=2.5 -> Non-integer subsampling factor `2.5`. A SubGrid's points must coincide with its parent's, so the factor must be a whole number of parent cells -- use `2` or `3`. Subsampling onto points that fall between parent points requires interpolation, which SubGrid does not perform\n", - "factor=0 -> Subsampling factor must be >= 1, got `0`\n", + "factor=0 -> Subsampling factor must be > 0, got `0`\n", + "factor=-1 -> Subsampling factor must be > 0, got `-1`\n", "factor=(2, 2, 2) -> Expected 1 or 2 factors for a 2D parent Grid, got 3\n" ] } ], "source": [ - "for bad in (2.5, 0, (2, 2, 2)):\n", + "for bad in (0, -1, (2, 2, 2)):\n", " try:\n", " SubGrid(grid, factor=bad)\n", " except ValueError as e:\n", @@ -365,10 +366,10 @@ "execution_count": 10, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.887762Z", - "iopub.status.busy": "2026-10-21T00:00:10.887135Z", - "iopub.status.idle": "2026-10-21T00:00:10.895555Z", - "shell.execute_reply": "2026-10-21T00:00:10.893836Z" + "iopub.execute_input": "2026-10-21T15:11:19.665029Z", + "iopub.status.busy": "2026-10-21T15:11:19.664542Z", + "iopub.status.idle": "2026-10-21T15:11:19.672124Z", + "shell.execute_reply": "2026-10-21T15:11:19.670540Z" } }, "outputs": [ @@ -402,10 +403,10 @@ "execution_count": 11, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.901337Z", - "iopub.status.busy": "2026-10-21T00:00:10.900718Z", - "iopub.status.idle": "2026-10-21T00:00:10.919067Z", - "shell.execute_reply": "2026-10-21T00:00:10.917441Z" + "iopub.execute_input": "2026-10-21T15:11:19.677450Z", + "iopub.status.busy": "2026-10-21T15:11:19.676871Z", + "iopub.status.idle": "2026-10-21T15:11:19.694697Z", + "shell.execute_reply": "2026-10-21T15:11:19.693284Z" } }, "outputs": [ @@ -439,10 +440,10 @@ "execution_count": 12, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.924379Z", - "iopub.status.busy": "2026-10-21T00:00:10.923847Z", - "iopub.status.idle": "2026-10-21T00:00:10.935092Z", - "shell.execute_reply": "2026-10-21T00:00:10.933488Z" + "iopub.execute_input": "2026-10-21T15:11:19.699405Z", + "iopub.status.busy": "2026-10-21T15:11:19.698986Z", + "iopub.status.idle": "2026-10-21T15:11:19.709783Z", + "shell.execute_reply": "2026-10-21T15:11:19.708238Z" } }, "outputs": [ @@ -478,10 +479,10 @@ "execution_count": 13, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.940307Z", - "iopub.status.busy": "2026-10-21T00:00:10.939772Z", - "iopub.status.idle": "2026-10-21T00:00:10.954034Z", - "shell.execute_reply": "2026-10-21T00:00:10.952426Z" + "iopub.execute_input": "2026-10-21T15:11:19.714884Z", + "iopub.status.busy": "2026-10-21T15:11:19.714394Z", + "iopub.status.idle": "2026-10-21T15:11:19.726918Z", + "shell.execute_reply": "2026-10-21T15:11:19.725499Z" } }, "outputs": [ @@ -521,10 +522,10 @@ "execution_count": 14, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:10.959450Z", - "iopub.status.busy": "2026-10-21T00:00:10.958997Z", - "iopub.status.idle": "2026-10-21T00:00:11.166399Z", - "shell.execute_reply": "2026-10-21T00:00:11.165173Z" + "iopub.execute_input": "2026-10-21T15:11:19.732185Z", + "iopub.status.busy": "2026-10-21T15:11:19.731617Z", + "iopub.status.idle": "2026-10-21T15:11:19.922943Z", + "shell.execute_reply": "2026-10-21T15:11:19.922123Z" } }, "outputs": [ @@ -575,10 +576,10 @@ "execution_count": 15, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.171003Z", - "iopub.status.busy": "2026-10-21T00:00:11.170781Z", - "iopub.status.idle": "2026-10-21T00:00:11.301506Z", - "shell.execute_reply": "2026-10-21T00:00:11.300355Z" + "iopub.execute_input": "2026-10-21T15:11:19.926929Z", + "iopub.status.busy": "2026-10-21T15:11:19.926760Z", + "iopub.status.idle": "2026-10-21T15:11:20.105763Z", + "shell.execute_reply": "2026-10-21T15:11:20.104394Z" } }, "outputs": [ @@ -586,6 +587,42 @@ "name": "stderr", "output_type": "stream", "text": [ + "Shrinking bounds (`x_m`, `x_M`); some `x` points will not be computed. Likely insufficient space_order for the derivatives.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:719: DevitoWarning: Shrinking bounds (`x_m`, `x_M`); some `x` points will not be computed. Likely insufficient space_order for the derivatives.\n", + " warn(f\"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); \"\n", + "Shrinking bounds (`xs_m`, `xs_M`); some `xs` points will not be computed. Likely insufficient space_order for the derivatives.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:719: DevitoWarning: Shrinking bounds (`xs_m`, `xs_M`); some `xs` points will not be computed. Likely insufficient space_order for the derivatives.\n", + " warn(f\"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); \"\n", + "Shrinking bounds (`y_m`, `y_M`); some `y` points will not be computed. Likely insufficient space_order for the derivatives.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:719: DevitoWarning: Shrinking bounds (`y_m`, `y_M`); some `y` points will not be computed. Likely insufficient space_order for the derivatives.\n", + " warn(f\"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); \"\n", + "Shrinking bounds (`ys_m`, `ys_M`); some `ys` points will not be computed. Likely insufficient space_order for the derivatives.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:719: DevitoWarning: Shrinking bounds (`ys_m`, `ys_M`); some `ys` points will not be computed. Likely insufficient space_order for the derivatives.\n", + " warn(f\"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); \"\n", "Operator `Kernel` ran in 0.01 s\n" ] }, @@ -629,10 +666,10 @@ "execution_count": 16, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.305798Z", - "iopub.status.busy": "2026-10-21T00:00:11.305576Z", - "iopub.status.idle": "2026-10-21T00:00:11.420274Z", - "shell.execute_reply": "2026-10-21T00:00:11.419116Z" + "iopub.execute_input": "2026-10-21T15:11:20.110313Z", + "iopub.status.busy": "2026-10-21T15:11:20.110034Z", + "iopub.status.idle": "2026-10-21T15:11:20.305359Z", + "shell.execute_reply": "2026-10-21T15:11:20.303927Z" } }, "outputs": [ @@ -676,10 +713,10 @@ "execution_count": 17, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.424454Z", - "iopub.status.busy": "2026-10-21T00:00:11.424226Z", - "iopub.status.idle": "2026-10-21T00:00:11.614968Z", - "shell.execute_reply": "2026-10-21T00:00:11.606198Z" + "iopub.execute_input": "2026-10-21T15:11:20.310118Z", + "iopub.status.busy": "2026-10-21T15:11:20.309765Z", + "iopub.status.idle": "2026-10-21T15:11:20.665081Z", + "shell.execute_reply": "2026-10-21T15:11:20.663606Z" } }, "outputs": [ @@ -731,10 +768,10 @@ "execution_count": 18, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.621540Z", - "iopub.status.busy": "2026-10-21T00:00:11.620971Z", - "iopub.status.idle": "2026-10-21T00:00:11.634243Z", - "shell.execute_reply": "2026-10-21T00:00:11.632756Z" + "iopub.execute_input": "2026-10-21T15:11:20.669884Z", + "iopub.status.busy": "2026-10-21T15:11:20.669587Z", + "iopub.status.idle": "2026-10-21T15:11:20.678615Z", + "shell.execute_reply": "2026-10-21T15:11:20.677321Z" } }, "outputs": [ @@ -765,10 +802,10 @@ "execution_count": 19, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.639205Z", - "iopub.status.busy": "2026-10-21T00:00:11.638777Z", - "iopub.status.idle": "2026-10-21T00:00:11.669891Z", - "shell.execute_reply": "2026-10-21T00:00:11.668737Z" + "iopub.execute_input": "2026-10-21T15:11:20.682805Z", + "iopub.status.busy": "2026-10-21T15:11:20.682561Z", + "iopub.status.idle": "2026-10-21T15:11:20.710333Z", + "shell.execute_reply": "2026-10-21T15:11:20.708992Z" } }, "outputs": [ @@ -803,10 +840,10 @@ "execution_count": 20, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.674050Z", - "iopub.status.busy": "2026-10-21T00:00:11.673852Z", - "iopub.status.idle": "2026-10-21T00:00:11.853934Z", - "shell.execute_reply": "2026-10-21T00:00:11.852728Z" + "iopub.execute_input": "2026-10-21T15:11:20.714735Z", + "iopub.status.busy": "2026-10-21T15:11:20.714462Z", + "iopub.status.idle": "2026-10-21T15:11:21.002353Z", + "shell.execute_reply": "2026-10-21T15:11:21.001018Z" } }, "outputs": [ @@ -864,10 +901,10 @@ "execution_count": 21, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:11.858435Z", - "iopub.status.busy": "2026-10-21T00:00:11.858182Z", - "iopub.status.idle": "2026-10-21T00:00:12.058665Z", - "shell.execute_reply": "2026-10-21T00:00:12.057133Z" + "iopub.execute_input": "2026-10-21T15:11:21.006682Z", + "iopub.status.busy": "2026-10-21T15:11:21.006432Z", + "iopub.status.idle": "2026-10-21T15:11:21.338580Z", + "shell.execute_reply": "2026-10-21T15:11:21.337310Z" } }, "outputs": [ @@ -917,10 +954,10 @@ "execution_count": 22, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.064134Z", - "iopub.status.busy": "2026-10-21T00:00:12.063653Z", - "iopub.status.idle": "2026-10-21T00:00:12.082047Z", - "shell.execute_reply": "2026-10-21T00:00:12.080721Z" + "iopub.execute_input": "2026-10-21T15:11:21.342886Z", + "iopub.status.busy": "2026-10-21T15:11:21.342622Z", + "iopub.status.idle": "2026-10-21T15:11:21.354554Z", + "shell.execute_reply": "2026-10-21T15:11:21.353299Z" } }, "outputs": [ @@ -945,8 +982,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Generic subsampling -- subgrid points that fall *between* parent points, requiring\n", - "interpolation -- is a separate topic and is not covered here." + "Generic subsampling -- subgrid points that fall *between* parent points -- **is**\n", + "covered, in the next section." ] }, { @@ -968,10 +1005,10 @@ "execution_count": 23, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.087228Z", - "iopub.status.busy": "2026-10-21T00:00:12.086738Z", - "iopub.status.idle": "2026-10-21T00:00:12.095143Z", - "shell.execute_reply": "2026-10-21T00:00:12.093635Z" + "iopub.execute_input": "2026-10-21T15:11:21.359212Z", + "iopub.status.busy": "2026-10-21T15:11:21.358948Z", + "iopub.status.idle": "2026-10-21T15:11:21.365140Z", + "shell.execute_reply": "2026-10-21T15:11:21.363935Z" } }, "outputs": [], @@ -998,10 +1035,10 @@ "execution_count": 24, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.100120Z", - "iopub.status.busy": "2026-10-21T00:00:12.099707Z", - "iopub.status.idle": "2026-10-21T00:00:12.239134Z", - "shell.execute_reply": "2026-10-21T00:00:12.237862Z" + "iopub.execute_input": "2026-10-21T15:11:21.369101Z", + "iopub.status.busy": "2026-10-21T15:11:21.368855Z", + "iopub.status.idle": "2026-10-21T15:11:21.569382Z", + "shell.execute_reply": "2026-10-21T15:11:21.567992Z" } }, "outputs": [ @@ -1051,10 +1088,10 @@ "execution_count": 25, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.244318Z", - "iopub.status.busy": "2026-10-21T00:00:12.243861Z", - "iopub.status.idle": "2026-10-21T00:00:12.252593Z", - "shell.execute_reply": "2026-10-21T00:00:12.251182Z" + "iopub.execute_input": "2026-10-21T15:11:21.573858Z", + "iopub.status.busy": "2026-10-21T15:11:21.573583Z", + "iopub.status.idle": "2026-10-21T15:11:21.580827Z", + "shell.execute_reply": "2026-10-21T15:11:21.579548Z" } }, "outputs": [ @@ -1089,10 +1126,10 @@ "execution_count": 26, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.258305Z", - "iopub.status.busy": "2026-10-21T00:00:12.257700Z", - "iopub.status.idle": "2026-10-21T00:00:12.384155Z", - "shell.execute_reply": "2026-10-21T00:00:12.382996Z" + "iopub.execute_input": "2026-10-21T15:11:21.585390Z", + "iopub.status.busy": "2026-10-21T15:11:21.585115Z", + "iopub.status.idle": "2026-10-21T15:11:21.777160Z", + "shell.execute_reply": "2026-10-21T15:11:21.775234Z" } }, "outputs": [ @@ -1127,10 +1164,10 @@ "execution_count": 27, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.388152Z", - "iopub.status.busy": "2026-10-21T00:00:12.387972Z", - "iopub.status.idle": "2026-10-21T00:00:12.393955Z", - "shell.execute_reply": "2026-10-21T00:00:12.392797Z" + "iopub.execute_input": "2026-10-21T15:11:21.782049Z", + "iopub.status.busy": "2026-10-21T15:11:21.781739Z", + "iopub.status.idle": "2026-10-21T15:11:21.790650Z", + "shell.execute_reply": "2026-10-21T15:11:21.789381Z" } }, "outputs": [ @@ -1161,10 +1198,10 @@ "execution_count": 28, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.398049Z", - "iopub.status.busy": "2026-10-21T00:00:12.397832Z", - "iopub.status.idle": "2026-10-21T00:00:12.520522Z", - "shell.execute_reply": "2026-10-21T00:00:12.518923Z" + "iopub.execute_input": "2026-10-21T15:11:21.795106Z", + "iopub.status.busy": "2026-10-21T15:11:21.794838Z", + "iopub.status.idle": "2026-10-21T15:11:22.171743Z", + "shell.execute_reply": "2026-10-21T15:11:22.170428Z" } }, "outputs": [ @@ -1211,10 +1248,10 @@ "execution_count": 29, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.525875Z", - "iopub.status.busy": "2026-10-21T00:00:12.525357Z", - "iopub.status.idle": "2026-10-21T00:00:12.755021Z", - "shell.execute_reply": "2026-10-21T00:00:12.753816Z" + "iopub.execute_input": "2026-10-21T15:11:22.176285Z", + "iopub.status.busy": "2026-10-21T15:11:22.176015Z", + "iopub.status.idle": "2026-10-21T15:11:22.376760Z", + "shell.execute_reply": "2026-10-21T15:11:22.375393Z" } }, "outputs": [ @@ -1249,10 +1286,10 @@ "execution_count": 30, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.758836Z", - "iopub.status.busy": "2026-10-21T00:00:12.758655Z", - "iopub.status.idle": "2026-10-21T00:00:12.889574Z", - "shell.execute_reply": "2026-10-21T00:00:12.887973Z" + "iopub.execute_input": "2026-10-21T15:11:22.381116Z", + "iopub.status.busy": "2026-10-21T15:11:22.380859Z", + "iopub.status.idle": "2026-10-21T15:11:22.573341Z", + "shell.execute_reply": "2026-10-21T15:11:22.572185Z" } }, "outputs": [ @@ -1309,10 +1346,10 @@ "execution_count": 31, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:12.894890Z", - "iopub.status.busy": "2026-10-21T00:00:12.894468Z", - "iopub.status.idle": "2026-10-21T00:00:13.027969Z", - "shell.execute_reply": "2026-10-21T00:00:13.026750Z" + "iopub.execute_input": "2026-10-21T15:11:22.577926Z", + "iopub.status.busy": "2026-10-21T15:11:22.577679Z", + "iopub.status.idle": "2026-10-21T15:11:22.766598Z", + "shell.execute_reply": "2026-10-21T15:11:22.765232Z" } }, "outputs": [ @@ -1350,10 +1387,10 @@ "execution_count": 32, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.031759Z", - "iopub.status.busy": "2026-10-21T00:00:13.031587Z", - "iopub.status.idle": "2026-10-21T00:00:13.051619Z", - "shell.execute_reply": "2026-10-21T00:00:13.050490Z" + "iopub.execute_input": "2026-10-21T15:11:22.770966Z", + "iopub.status.busy": "2026-10-21T15:11:22.770679Z", + "iopub.status.idle": "2026-10-21T15:11:22.795688Z", + "shell.execute_reply": "2026-10-21T15:11:22.794390Z" } }, "outputs": [ @@ -1410,10 +1447,10 @@ "execution_count": 33, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.055748Z", - "iopub.status.busy": "2026-10-21T00:00:13.055580Z", - "iopub.status.idle": "2026-10-21T00:00:13.189811Z", - "shell.execute_reply": "2026-10-21T00:00:13.188606Z" + "iopub.execute_input": "2026-10-21T15:11:22.800339Z", + "iopub.status.busy": "2026-10-21T15:11:22.800075Z", + "iopub.status.idle": "2026-10-21T15:11:23.020059Z", + "shell.execute_reply": "2026-10-21T15:11:23.018704Z" } }, "outputs": [ @@ -1454,10 +1491,10 @@ "execution_count": 34, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.193859Z", - "iopub.status.busy": "2026-10-21T00:00:13.193626Z", - "iopub.status.idle": "2026-10-21T00:00:13.218664Z", - "shell.execute_reply": "2026-10-21T00:00:13.217527Z" + "iopub.execute_input": "2026-10-21T15:11:23.024463Z", + "iopub.status.busy": "2026-10-21T15:11:23.024170Z", + "iopub.status.idle": "2026-10-21T15:11:23.055530Z", + "shell.execute_reply": "2026-10-21T15:11:23.054230Z" } }, "outputs": [ @@ -1472,8 +1509,6 @@ "name": "stderr", "output_type": "stream", "text": [ - "/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:718: DevitoWarning: Shrinking bounds (`x_m`, `x_M`); some `x` points will not be computed. Likely insufficient space_order for the derivatives.\n", - " warn(f\"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); \"\n", "Shrinking bounds (`xs_m`, `xs_M`); some `xs` points will not be computed. Likely insufficient space_order for the derivatives.\n" ] }, @@ -1481,8 +1516,6 @@ "name": "stderr", "output_type": "stream", "text": [ - "/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:718: DevitoWarning: Shrinking bounds (`xs_m`, `xs_M`); some `xs` points will not be computed. Likely insufficient space_order for the derivatives.\n", - " warn(f\"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); \"\n", "Operator `Kernel` ran in 0.01 s\n" ] }, @@ -1522,10 +1555,10 @@ "execution_count": 35, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.222737Z", - "iopub.status.busy": "2026-10-21T00:00:13.222568Z", - "iopub.status.idle": "2026-10-21T00:00:13.400826Z", - "shell.execute_reply": "2026-10-21T00:00:13.399620Z" + "iopub.execute_input": "2026-10-21T15:11:23.060091Z", + "iopub.status.busy": "2026-10-21T15:11:23.059805Z", + "iopub.status.idle": "2026-10-21T15:11:23.364064Z", + "shell.execute_reply": "2026-10-21T15:11:23.362647Z" } }, "outputs": [ @@ -1590,10 +1623,10 @@ "execution_count": 36, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.404763Z", - "iopub.status.busy": "2026-10-21T00:00:13.404595Z", - "iopub.status.idle": "2026-10-21T00:00:13.428563Z", - "shell.execute_reply": "2026-10-21T00:00:13.427431Z" + "iopub.execute_input": "2026-10-21T15:11:23.368455Z", + "iopub.status.busy": "2026-10-21T15:11:23.368162Z", + "iopub.status.idle": "2026-10-21T15:11:23.400374Z", + "shell.execute_reply": "2026-10-21T15:11:23.399047Z" } }, "outputs": [ @@ -1641,10 +1674,10 @@ "execution_count": 37, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.432786Z", - "iopub.status.busy": "2026-10-21T00:00:13.432609Z", - "iopub.status.idle": "2026-10-21T00:00:13.568839Z", - "shell.execute_reply": "2026-10-21T00:00:13.567306Z" + "iopub.execute_input": "2026-10-21T15:11:23.404920Z", + "iopub.status.busy": "2026-10-21T15:11:23.404665Z", + "iopub.status.idle": "2026-10-21T15:11:23.641038Z", + "shell.execute_reply": "2026-10-21T15:11:23.639634Z" } }, "outputs": [ @@ -1705,10 +1738,10 @@ "execution_count": 38, "metadata": { "execution": { - "iopub.execute_input": "2026-10-21T00:00:13.574150Z", - "iopub.status.busy": "2026-10-21T00:00:13.573710Z", - "iopub.status.idle": "2026-10-21T00:00:13.750574Z", - "shell.execute_reply": "2026-10-21T00:00:13.749383Z" + "iopub.execute_input": "2026-10-21T15:11:23.645768Z", + "iopub.status.busy": "2026-10-21T15:11:23.645496Z", + "iopub.status.idle": "2026-10-21T15:11:23.935804Z", + "shell.execute_reply": "2026-10-21T15:11:23.934381Z" } }, "outputs": [ @@ -1757,39 +1790,609 @@ "show(op_indirect)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interpolated subgrids\n", + "\n", + "A non-integer factor puts subgrid points *between* parent points, so reading such a\n", + "Function on the parent grid has to interpolate. This needs no new API: a SubGrid\n", + "Function used in an equation is evaluated **at the left-hand side's location**, exactly\n", + "as a staggered Function is." + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:23.940688Z", + "iopub.status.busy": "2026-10-21T15:11:23.940420Z", + "iopub.status.idle": "2026-10-21T15:11:23.950991Z", + "shell.execute_reply": "2026-10-21T15:11:23.949727Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "parent shape (13, 13) spacing (1.0, 1.0)\n", + "subgrid shape (9, 9) spacing (1.5, 1.5)\n", + "interpolated True\n", + "dimensions (xs, ys)\n" + ] + } + ], + "source": [ + "grid_i = Grid(shape=(13, 13), extent=(12., 12.))\n", + "sub_i = SubGrid(grid_i, factor=1.5)\n", + "\n", + "spacing = lambda g: tuple(float(s) for s in g.spacing)\n", + "print('parent shape ', grid_i.shape, ' spacing', spacing(grid_i))\n", + "print('subgrid shape', sub_i.shape, ' spacing', spacing(sub_i))\n", + "print('interpolated ', sub_i.is_interpolated)\n", + "print('dimensions ', sub_i.dimensions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`dimensions` are now plain `SpaceDimension`s with their own spacing, not\n", + "`ConditionalDimension`s: there is no stride to exploit, because the points do not line\n", + "up with the parent's.\n", + "\n", + "### Reading a subgrid Function on the parent\n", + "\n", + "`Eq(u, vp)` reads `vp` at `u`'s points -- no wrapper and no flag. The relocation is\n", + "recorded symbolically at first and expanded into a stencil during compilation, so the\n", + "interpolation weights can be chosen once the iteration space is known." + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:23.955139Z", + "iopub.status.busy": "2026-10-21T15:11:23.954897Z", + "iopub.status.idle": "2026-10-21T15:11:23.966265Z", + "shell.execute_reply": "2026-10-21T15:11:23.965006Z" + } + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle u_i(x, y) = \\operatorname{SubGridInterp}\\left(vp_i(xs, ys), u_i(x, y)\\right)$" + ], + "text/plain": [ + "Eq(u_i(x, y), SubGridInterp(vp_i(xs, ys), at=u_i(x, y)))" + ] + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u_i = Function(name='u_i', grid=grid_i, space_order=2)\n", + "vp_i = Function(name='vp_i', grid=sub_i, space_order=2)\n", + "\n", + "eq = Eq(u_i, vp_i)\n", + "eq.evaluate" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`SubGridInterp(vp_i(xs, ys), at=u_i(x, y))` is that record: *read `vp_i` at `u_i`'s\n", + "location*. The compiler turns it into the stencil below.\n", + "\n", + "### The generated code\n", + "\n", + "The weights depend on the parent index only through `x % 3`, so they live in a small\n", + "compile-time table and the loops stay plain and unit-stride -- no branch, no `floorf`,\n", + "and only the parent's Dimensions are iterated." + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:23.970437Z", + "iopub.status.busy": "2026-10-21T15:11:23.970178Z", + "iopub.status.idle": "2026-10-21T15:11:24.218910Z", + "shell.execute_reply": "2026-10-21T15:11:24.217568Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Devito generated code for Operator `Kernel` */\n", + "\n", + "#define _POSIX_C_SOURCE 200809L\n", + "#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);\n", + "#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;\n", + "\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct dataobj\n", + "{\n", + " void *restrict data;\n", + " int * size;\n", + " unsigned long nbytes;\n", + " unsigned long * npsize;\n", + " unsigned long * dsize;\n", + " int * hsize;\n", + " int * hofs;\n", + " int * oofs;\n", + " void * dmap;\n", + "} ;\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(struct dataobj *restrict u_i_vec, struct dataobj *restrict vp_i_vec, const int x_M, const int x_m, const int y_M, const int y_m, struct profiler * timers)\n", + "{\n", + " float w0[6] __attribute__ ((aligned (64))) = {1.0F, 0.0F, 3.33333333333333e-1F, 6.66666666666667e-1F, 6.66666666666667e-1F, 3.33333333333333e-1F};\n", + "\n", + " float (*restrict u_i)[u_i_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[u_i_vec->size[1]]) u_i_vec->data;\n", + " float (*restrict vp_i)[vp_i_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[vp_i_vec->size[1]]) vp_i_vec->data;\n", + "\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + "\n", + " START(section0)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " #pragma omp simd aligned(u_i,vp_i:64)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " u_i[x + 2][y + 2] = w0[2*((x)%(3)) + 1]*w0[2*((y)%(3)) + 1]*vp_i[((2*x) / 3) + 3][((2*y) / 3) + 3] + w0[2*((x)%(3)) + 1]*w0[2*((y)%(3))]*vp_i[((2*x) / 3) + 3][((2*y) / 3) + 2] + w0[2*((y)%(3)) + 1]*w0[2*((x)%(3))]*vp_i[((2*x) / 3) + 2][((2*y) / 3) + 3] + w0[2*((x)%(3))]*w0[2*((y)%(3))]*vp_i[((2*x) / 3) + 2][((2*y) / 3) + 2];\n", + " }\n", + " }\n", + " STOP(section0,timers)\n", + "\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "op_i = Operator(eq)\n", + "print(op_i)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Reading at parent point `x` evaluates at subgrid coordinate `x/factor`. For\n", + "`factor = p/q` that is `x*q/p`, so the base index is `(x*q)//p` and the fractional part\n", + "`((x*q) % p)/p` repeats with period `p` -- hence a table of `p` phases, two linear taps\n", + "each.\n", + "\n", + "### Linear interpolation is exact on a linear field\n", + "\n", + "The sharpest check available: any error in a base index or a weight shows up at once." + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:24.223511Z", + "iopub.status.busy": "2026-10-21T15:11:24.223256Z", + "iopub.status.idle": "2026-10-21T15:11:24.253166Z", + "shell.execute_reply": "2026-10-21T15:11:24.251894Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "max error: 3.178914390389309e-06\n" + ] + } + ], + "source": [ + "idx = np.meshgrid(*[np.arange(n) for n in sub_i.shape], indexing='ij')\n", + "vp_i.data[:] = (3.*idx[0] + 2.*idx[1]).astype(vp_i.dtype)\n", + "u_i.data[:] = 0.\n", + "\n", + "op_i.apply()\n", + "\n", + "pidx = np.meshgrid(*[np.arange(n) for n in grid_i.shape], indexing='ij')\n", + "exact = 3.*(pidx[0]/1.5) + 2.*(pidx[1]/1.5)\n", + "\n", + "# Only points whose two taps stay inside the subgrid\n", + "inside = (pidx[0] <= int((sub_i.shape[0] - 1)*1.5) - 1) & \\\n", + " (pidx[1] <= int((sub_i.shape[1] - 1)*1.5) - 1)\n", + "print('max error:', np.abs(np.asarray(u_i.data)[inside] - exact[inside]).max())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### A factor below 1 refines\n", + "\n", + "Nothing in the lowering assumes coarsening -- `factor=0.5` makes the subgrid *finer*\n", + "than its parent. Every parent point then lands exactly on a subgrid point, so there is\n", + "a single phase: no weight table is emitted at all and the read collapses to a direct\n", + "strided copy." + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:24.257489Z", + "iopub.status.busy": "2026-10-21T15:11:24.257242Z", + "iopub.status.idle": "2026-10-21T15:11:24.461955Z", + "shell.execute_reply": "2026-10-21T15:11:24.460597Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "subgrid shape (25, 25) vs parent (13, 13)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Devito generated code for Operator `Kernel` */\n", + "\n", + "#define _POSIX_C_SOURCE 200809L\n", + "#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);\n", + "#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;\n", + "\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct dataobj\n", + "{\n", + " void *restrict data;\n", + " int * size;\n", + " unsigned long nbytes;\n", + " unsigned long * npsize;\n", + " unsigned long * dsize;\n", + " int * hsize;\n", + " int * hofs;\n", + " int * oofs;\n", + " void * dmap;\n", + "} ;\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(struct dataobj *restrict u_i_vec, struct dataobj *restrict vp_f_vec, const int x_M, const int x_m, const int y_M, const int y_m, struct profiler * timers)\n", + "{\n", + " float (*restrict u_i)[u_i_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[u_i_vec->size[1]]) u_i_vec->data;\n", + " float (*restrict vp_f)[vp_f_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[vp_f_vec->size[1]]) vp_f_vec->data;\n", + "\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + "\n", + " START(section0)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " #pragma omp simd aligned(u_i,vp_f:64)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " u_i[x + 2][y + 2] = vp_f[2*x + 2][2*y + 2];\n", + " }\n", + " }\n", + " STOP(section0,timers)\n", + "\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "sub_f = SubGrid(grid_i, factor=0.5)\n", + "print('subgrid shape', sub_f.shape, 'vs parent', grid_i.shape)\n", + "\n", + "vp_f = Function(name='vp_f', grid=sub_f, space_order=2)\n", + "print(Operator(Eq(u_i, vp_f)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Per-axis factors\n", + "\n", + "Each axis carries its own factor. A factor of exactly `1` passes the parent's own\n", + "Dimension straight through, so that axis is neither strided nor interpolated.\n", + "\n", + "Note that *within* an interpolated SubGrid even an integer axis interpolates. With\n", + "`factor=2`, the parent point `x=1` lies halfway between two subgrid points, so reading\n", + "it genuinely needs both taps -- unlike the pure subsampling case, where the subgrid is\n", + "written and only its own points are visited." + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:24.466465Z", + "iopub.status.busy": "2026-10-21T15:11:24.466170Z", + "iopub.status.idle": "2026-10-21T15:11:24.749143Z", + "shell.execute_reply": "2026-10-21T15:11:24.747773Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dimensions (xs, ys, z) shape (7, 9, 13)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Devito generated code for Operator `Kernel` */\n", + "\n", + "#define _POSIX_C_SOURCE 200809L\n", + "#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);\n", + "#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;\n", + "\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct dataobj\n", + "{\n", + " void *restrict data;\n", + " int * size;\n", + " unsigned long nbytes;\n", + " unsigned long * npsize;\n", + " unsigned long * dsize;\n", + " int * hsize;\n", + " int * hofs;\n", + " int * oofs;\n", + " void * dmap;\n", + "} ;\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(struct dataobj *restrict u_m_vec, struct dataobj *restrict vp_m_vec, const int x_M, const int x_m, const int y_M, const int y_m, const int z_M, const int z_m, struct profiler * timers)\n", + "{\n", + " float w0[4] __attribute__ ((aligned (64))) = {1.0F, 0.0F, 5.0e-1F, 5.0e-1F};\n", + " float w1[6] __attribute__ ((aligned (64))) = {1.0F, 0.0F, 3.33333333333333e-1F, 6.66666666666667e-1F, 6.66666666666667e-1F, 3.33333333333333e-1F};\n", + "\n", + " float (*restrict u_m)[u_m_vec->size[1]][u_m_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[u_m_vec->size[1]][u_m_vec->size[2]]) u_m_vec->data;\n", + " float (*restrict vp_m)[vp_m_vec->size[1]][vp_m_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[vp_m_vec->size[1]][vp_m_vec->size[2]]) vp_m_vec->data;\n", + "\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + "\n", + " START(section0)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " #pragma omp simd aligned(u_m,vp_m:64)\n", + " for (int z = z_m; z <= z_M; z += 1)\n", + " {\n", + " u_m[x + 2][y + 2][z + 2] = w0[2*((x & 1)) + 1]*w1[2*((y)%(3)) + 1]*vp_m[(x / 2) + 3][((2*y) / 3) + 3][z + 2] + w0[2*((x & 1)) + 1]*w1[2*((y)%(3))]*vp_m[(x / 2) + 3][((2*y) / 3) + 2][z + 2] + w0[2*((x & 1))]*w1[2*((y)%(3)) + 1]*vp_m[(x / 2) + 2][((2*y) / 3) + 3][z + 2] + w0[2*((x & 1))]*w1[2*((y)%(3))]*vp_m[(x / 2) + 2][((2*y) / 3) + 2][z + 2];\n", + " }\n", + " }\n", + " }\n", + " STOP(section0,timers)\n", + "\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "grid_m = Grid(shape=(13, 13, 13), extent=(12.,)*3)\n", + "sub_m = SubGrid(grid_m, factor=(2, 1.5, 1.0))\n", + "print('dimensions', sub_m.dimensions, ' shape', sub_m.shape)\n", + "\n", + "u_m = Function(name='u_m', grid=grid_m, space_order=2)\n", + "vp_m = Function(name='vp_m', grid=sub_m, space_order=2)\n", + "print(Operator(Eq(u_m, vp_m)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One weight table per axis that needs interpolation (`w0` for `x`, 2 phases; `w1` for\n", + "`y`, 3 phases) and none for the `z` axis of factor 1. All three loops are plain and\n", + "unit-stride.\n", + "\n", + "### Integer-factor SubGrids still subsample\n", + "\n", + "A SubGrid whose factors are *all* integer keeps the Part 1 behaviour -- a strided loop,\n", + "no interpolation and no weights -- so the two regimes coexist." + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": { + "execution": { + "iopub.execute_input": "2026-10-21T15:11:24.753796Z", + "iopub.status.busy": "2026-10-21T15:11:24.753457Z", + "iopub.status.idle": "2026-10-21T15:11:24.966668Z", + "shell.execute_reply": "2026-10-21T15:11:24.965316Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "interpolated: False dimensions: (xk, yk)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Devito generated code for Operator `Kernel` */\n", + "\n", + "#define _POSIX_C_SOURCE 200809L\n", + "#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);\n", + "#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;\n", + "\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct dataobj\n", + "{\n", + " void *restrict data;\n", + " int * size;\n", + " unsigned long nbytes;\n", + " unsigned long * npsize;\n", + " unsigned long * dsize;\n", + " int * hsize;\n", + " int * hofs;\n", + " int * oofs;\n", + " void * dmap;\n", + "} ;\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(struct dataobj *restrict u_i_vec, struct dataobj *restrict vp_int_vec, const int xk_M, const int xk_m, const int xkf, const int yk_M, const int yk_m, const int ykf, struct profiler * timers)\n", + "{\n", + " float (*restrict u_i)[u_i_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[u_i_vec->size[1]]) u_i_vec->data;\n", + " float (*restrict vp_int)[vp_int_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[vp_int_vec->size[1]]) vp_int_vec->data;\n", + "\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + "\n", + " START(section0)\n", + " for (int xk = xk_m; xk <= xk_M; xk += 1)\n", + " {\n", + " #pragma omp simd aligned(u_i,vp_int:64)\n", + " for (int yk = yk_m; yk <= yk_M; yk += 1)\n", + " {\n", + " vp_int[xk + 2][yk + 2] = u_i[xk*xkf + 2][yk*ykf + 2];\n", + " }\n", + " }\n", + " STOP(section0,timers)\n", + "\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "sub_int = SubGrid(grid_i, factor=3, name='k')\n", + "print('interpolated:', sub_int.is_interpolated, ' dimensions:', sub_int.dimensions)\n", + "\n", + "vp_int = Function(name='vp_int', grid=sub_int, space_order=2)\n", + "print(Operator(Eq(vp_int, u_i)))" + ] + }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", - "`SubGrid(parent, factor)` builds a coarser grid from a parent `Grid`, keeping every\n", - "`factor`-th point. It is a `Grid` subclass, so Functions, TimeFunctions, Operators and\n", - "finite differences all work as usual:\n", + "`SubGrid(parent, factor)` builds a grid from a parent `Grid`. It is a `Grid` subclass,\n", + "so Functions, TimeFunctions, Operators and finite differences all work as usual, and\n", + "the factor decides which of two regimes applies.\n", "\n", "```python\n", - "sub = SubGrid(grid, factor=2) # or factor=(2, 3); 1 = passthrough\n", - "usub = Function(name='usub', grid=sub) # shape and dimensions derived\n", + "sub = SubGrid(grid, factor=2) # integer -> subsampling\n", + "usub = Function(name='usub', grid=sub)\n", "Operator(Eq(usub, u)).apply() # strided copy from the parent\n", + "\n", + "sub = SubGrid(grid, factor=1.5) # non-integer -> interpolation\n", + "vp = Function(name='vp', grid=sub)\n", + "Operator(Eq(u, vp*u)).apply() # vp read at u's points\n", "```\n", "\n", - "| | derived from the parent |\n", - "|---|---|\n", - "| `shape` | `ceil(n/factor)` per axis |\n", - "| `spacing` | `factor * parent.spacing` |\n", - "| `extent`, `origin`, `dtype`, `time_dim` | inherited unchanged |\n", - "| `dimensions` | `ConditionalDimension`s on the parent's |\n", - "| `parent`, `factor` | the two attributes `SubGrid` adds |\n", + "| | integer factor | non-integer factor |\n", + "|---|---|---|\n", + "| `is_interpolated` | `False` | `True` |\n", + "| `shape` | `ceil(n/factor)` | `floor((n-1)/factor) + 1` |\n", + "| `dimensions` | `ConditionalDimension`s on the parent's | own `SpaceDimension`s |\n", + "| reading on the parent | strided, exact | interpolated at the LHS |\n", + "\n", + "Per-axis factors are independent, and `1` passes the parent's own Dimension through\n", + "unchanged. Common to both: `spacing` is `factor * parent.spacing`; `extent`, `origin`,\n", + "`dtype` and `time_dim` are inherited; `parent` and `factor` are the attributes\n", + "`SubGrid` adds.\n", "\n", "Use a `ConditionalDimension` directly for subsampling in **time**, or to combine a\n", "factor with an extra `condition` -- `SubGrid` covers neither.\n", "\n", - "Under the hood Devito iterates the subgrid with a unit stride and rescales parent\n", - "accesses, falling back to a pointwise `if (x%factor == 0)` guard in the four cases\n", - "above. Both compute the same thing; only the generated loop differs.\n", - "\n", - "Generic subsampling -- subgrid points that fall *between* parent grid points and so\n", - "require interpolation -- is a separate topic, not covered here." + "Under the hood, an integer factor iterates the subgrid with a unit stride and rescales\n", + "parent accesses, falling back to a pointwise `if (x%factor == 0)` guard in the four\n", + "cases above. A non-integer factor becomes a plain two-tap stencil whose weights come\n", + "from a compile-time table indexed by the phase `x % p`, so the loops stay unit-stride\n", + "and branch-free. Linear interpolation is used today; it is exact on a linear field.\n" ] } ], diff --git a/tests/test_subgrid.py b/tests/test_subgrid.py new file mode 100644 index 0000000000..51b00c91df --- /dev/null +++ b/tests/test_subgrid.py @@ -0,0 +1,288 @@ +import numpy as np +import pytest + +from devito import Eq, Function, Grid, Operator, SubGrid, switchconfig +from devito.ir.iet import Expression, FindNodes, Iteration +from devito.symbolics import retrieve_indexed + +COEFFS = (3., 2., 0.5) + + +def linear(shape): + """`3*i + 2*j + 0.5*k` over `shape`, in index units.""" + idx = np.meshgrid(*[np.arange(n) for n in shape], indexing='ij') + return sum(c*i for c, i in zip(COEFFS[:len(shape)], idx, strict=True)) + + +def restricted(factor, n=25, ndim=2): + """Read a linear SubGrid field at its parent's points.""" + grid = Grid(shape=(n,)*ndim, extent=(float(n - 1),)*ndim) + sub = SubGrid(grid, factor=factor) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + vp.data[:] = linear(sub.shape).astype(vp.dtype) + u.data[:] = 0. + + Operator(Eq(u, vp)).apply() + + return sub, u, vp + + +def interior(sub, n, ndim): + """Mask of parent points whose stencil stays inside the SubGrid.""" + idx = np.meshgrid(*[np.arange(n)]*ndim, indexing='ij') + mask = np.ones_like(idx[0], dtype=bool) + for i, (ns, f) in enumerate(zip(sub.shape, sub.factor, strict=True)): + mask &= idx[i] <= int((ns - 1)*float(f)) - 1 + return mask, idx + + +class TestApi: + + """ + A SubGrid Function is read at the LHS's location, with no wrapper -- the + same convention a staggered Function follows. + """ + + def test_plain_equation(self): + grid = Grid(shape=(13, 13), extent=(12., 12.)) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + u.data[:] = 1. + vp.data[:] = 2.5 + + Operator(Eq(u, vp*u)).apply() + + mask, _ = interior(sub, 13, 2) + assert np.allclose(np.asarray(u.data)[mask], 2.5, atol=1e-5) + + def test_iterates_parent_dimensions_only(self): + """ + The SubGrid Dimensions are not iterated: `vp` is read through a + stencil indexed by the parent's Dimensions. + """ + grid = Grid(shape=(13, 13), extent=(12., 12.)) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + + op = Operator(Eq(u, vp)) + + dims = [i.dim.name for i in FindNodes(Iteration).visit(op)] + assert dims == [d.name for d in grid.dimensions] + + def test_every_axis_relocated_under_blocking(self): + """ + In 3D, blocking has already run when the subgrid pass sees the Cluster, + so a parent Dimension appears in the IterationSpace only through its + block Dimensions. Every axis must still be relocated -- the symptom of + getting this wrong is a read left indexed by a SubGrid Dimension. + """ + grid = Grid(shape=(17, 17, 17), extent=(16.,)*3) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + + op = Operator(Eq(u, vp)) + + e, = [i.expr for i in FindNodes(Expression).visit(op) if i.write is u] + reads = [i for i in retrieve_indexed(e) if i.function is vp] + assert reads + + gone = set(sub.dimensions) + for i in reads: + assert not (set(i.free_symbols) & gone) + + # And the loop nest is the parent's, not a deeper one + dims = [i.dim for i in FindNodes(Iteration).visit(op)] + assert not (set(dims) & gone) + + def test_shared_factor_one_dimension_still_iterates(self): + """ + An axis of factor 1 shares the parent's own Dimension, so it must stay + in the IterationSpace -- dropping it leaves a loop bound with no + runtime value. + """ + grid = Grid(shape=(13, 13, 13), extent=(12.,)*3) + sub = SubGrid(grid, factor=(2, 1.5, 1.0)) + assert sub.dimensions[2] is grid.dimensions[2] + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + vp.data[:] = 2.5 + + op = Operator(Eq(u, vp)) + + dims = [i.dim.name for i in FindNodes(Iteration).visit(op)] + assert dims == [d.name for d in grid.dimensions] + + op.apply() # must not raise: `z` needs its bounds + assert np.allclose(np.asarray(u.data)[:-1, :-1, :], 2.5, atol=1e-5) + + def test_is_a_plain_stencil(self): + """ + No accumulator and no per-tap branch: the taps are a flat weighted + sum, so the inner loop stays vectorizable. + """ + grid = Grid(shape=(25, 25), extent=(24., 24.)) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + + with switchconfig(language='openmp'): + ccode = str(Operator(Eq(u, vp)).ccode) + + assert 'for (int r' not in ccode + assert 'if (' not in ccode + assert 'floor' not in ccode + assert '#pragma omp simd' in ccode + + def test_weights_are_compile_time(self): + """The per-phase weights are a stack table, not computed per point.""" + grid = Grid(shape=(13, 13), extent=(12., 12.)) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + + ccode = str(Operator(Eq(u, vp)).ccode) + + # `factor = 3/2` reading at the parent's points gives 3 phases of 2 + # taps, i.e. a 6-entry table + assert '[6]' in ccode + + @pytest.mark.parametrize('ndim,taps', [(1, 2), (2, 4), (3, 8)]) + def test_tap_count(self, ndim, taps): + """Linear interpolation reads `2**ndim` SubGrid points.""" + n = 13 + grid = Grid(shape=(n,)*ndim, extent=(float(n - 1),)*ndim) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + + op = Operator(Eq(u, vp)) + expr, = [i.expr for i in FindNodes(Expression).visit(op) if i.write is u] + + assert len([i for i in retrieve_indexed(expr.rhs) + if i.function is vp]) == taps + + +class TestExactness: + + """ + Linear interpolation is exact on a linear field -- the sharpest check, as + any error in a base index or a weight shows up immediately. + """ + + @pytest.mark.parametrize('factor', [1.5, 2.5, 0.5, 0.75, 1.25]) + def test_1d(self, factor): + n = 25 + sub, u, _ = restricted(factor, n=n, ndim=1) + + mask, idx = interior(sub, n, 1) + expected = 3.*(idx[0]/float(sub.factor[0])) + + assert np.allclose(np.asarray(u.data)[mask], expected[mask], atol=2e-4) + + @pytest.mark.parametrize('factor', [1.5, 2.5, 0.5, (2, 1.5), (0.5, 1.5)]) + def test_2d(self, factor): + n = 25 + sub, u, _ = restricted(factor, n=n, ndim=2) + + mask, idx = interior(sub, n, 2) + f0, f1 = [float(v) for v in sub.factor] + expected = 3.*(idx[0]/f0) + 2.*(idx[1]/f1) + + assert np.allclose(np.asarray(u.data)[mask], expected[mask], atol=2e-4) + + @pytest.mark.parametrize('factor', [1.5, 2.5, (2, 1.5, 1.)]) + def test_3d(self, factor): + """ + 3D is the case where blocking has already run by the time the subgrid + pass sees the Cluster, so a parent Dimension shows up only through its + block Dimensions. Every axis must still be relocated. + """ + n = 17 + sub, u, _ = restricted(factor, n=n, ndim=3) + + mask, idx = interior(sub, n, 3) + expected = sum(c*(i/float(f)) for c, i, f in + zip(COEFFS, idx, sub.factor, strict=True)) + + assert np.allclose(np.asarray(u.data)[mask], expected[mask], atol=2e-4) + + def test_constant_field(self): + """The weights must sum to one at every phase.""" + n = 25 + grid = Grid(shape=(n, n), extent=(float(n - 1),)*2) + sub = SubGrid(grid, factor=1.5) + + u = Function(name='u', grid=grid, space_order=2) + vp = Function(name='vp', grid=sub, space_order=2) + vp.data[:] = -3.25 + u.data[:] = 0. + + Operator(Eq(u, vp)).apply() + + mask, _ = interior(sub, n, 2) + assert np.allclose(np.asarray(u.data)[mask], -3.25, atol=1e-5) + + def test_refinement(self): + """A factor below 1 refines: every parent point is an exact hit.""" + n = 25 + sub, u, vp = restricted(0.5, n=n, ndim=1) + + assert sub.shape[0] > n + mask, idx = interior(sub, n, 1) + assert np.allclose(np.asarray(u.data)[mask], (3.*idx[0]/0.5)[mask], + atol=1e-4) + + +class TestIntegerFactor: + + """An integer factor needs no interpolation and keeps Part 1's lowering.""" + + def test_not_interpolated(self): + grid = Grid(shape=(12, 12)) + + assert not SubGrid(grid, factor=3).is_interpolated + assert SubGrid(grid, factor=1.5).is_interpolated + + def test_strided_copy(self): + n = 24 + grid = Grid(shape=(n, n)) + sub = SubGrid(grid, factor=3) + + u = Function(name='u', grid=grid) + g = Function(name='g', grid=sub, space_order=0) + data = np.arange(n*n, dtype=u.dtype).reshape(n, n) + u.data[:] = data + + Operator(Eq(g, u)).apply() + + assert np.array_equal(np.asarray(g.data), data[::3, ::3]) + + def test_self_reference(self): + """ + A SubGrid Function read at a shift of its own Dimensions is a + recurrence, not a relocation -- it must not be interpolated. + """ + grid = Grid(shape=(12, 12)) + sub = SubGrid(grid, factor=(3, 1)) + xs = sub.dimensions[0] + + f = Function(name='f', grid=grid) + g = Function(name='g', grid=sub, space_order=0) + f.data[:] = 1. + g.data[:] = 0. + + Operator(Eq(g, g.subs(xs, xs - 1) + f)).apply() + + assert np.allclose(np.asarray(g.data)[:, 0], [0., 1., 2., 3.]) diff --git a/tests/test_subsampling.py b/tests/test_subsampling.py index 07ff1a5a1f..7e1dcc006e 100644 --- a/tests/test_subsampling.py +++ b/tests/test_subsampling.py @@ -1,5 +1,3 @@ -import re - import numpy as np import pytest from sympy import Or @@ -734,9 +732,8 @@ def test_derivative_numerics(self): assert np.allclose(d.data[1:-1, 1:-1], 3., atol=1e-4) @pytest.mark.parametrize('factor,err', [ - (0, 'must be >= 1'), - (-1, 'must be >= 1'), - (2.5, 'Non-integer'), + (0, 'must be > 0'), + (-1, 'must be > 0'), ((2, 2, 2), 'Expected 1 or 2 factors'), ]) def test_invalid_factor(self, factor, err): @@ -791,18 +788,15 @@ def test_nesting_suggestion_is_equivalent(self, f0, f1, combined): assert np.array_equal(f.data, two_step) - def test_non_integer_factor_suggests_valid_factors(self): - """The suggested fallbacks must themselves be legal factors.""" + def test_non_integer_factor_is_interpolated(self): + """ + A non-integer factor is legal; it puts the SubGrid in the interpolated + regime, covered by `tests/test_subgrid.py`. + """ grid = Grid(shape=(9, 9)) - for bad in (0.5, 1.2, 2.5): - with pytest.raises(ValueError) as e: - SubGrid(grid, factor=bad) - suggested = re.findall(r'use `(\d+)` or `(\d+)`', str(e.value)) - assert suggested, str(e.value) - for v in suggested[0]: - assert int(v) >= 1 - SubGrid(grid, factor=int(v)) # must not raise + assert SubGrid(grid, factor=1.5).is_interpolated + assert not SubGrid(grid, factor=2).is_interpolated def test_parent_must_be_grid(self): with pytest.raises(ValueError, match='must be a Grid'):