compiler: Fix issue #2235

[georgebisbas]

In place of #2933.

The bug (issue #2235)

grid = Grid(shape=(4, 4))
u = TimeFunction(name='u', grid=grid)                 # regular, save=None
usave = TimeFunction(name='usave', grid=grid, save=5)  # save=5

op = Operator([Eq(u.forward, u + 1), Eq(usave, u)])
assert op.arguments()['time_M'] == 4   # actual: 3

time_M comes out one too small whenever a save-mode TimeFunction is
mixed with a regular (modulo-buffered) one.

Root cause

Cluster.dspace builds a per-Function IntervalGroup of access offsets
(parts), then unions them into one global Interval per Dimension:

intervals = IntervalGroup.generate('union', *parts.values())

Inspecting op._dspace.parts for the reproducer above:

u     : IntervalGroup[time[0,1], x[1,1], y[1,1]]   f.dimensions = (t, x, y)
usave : IntervalGroup[time[0,0], x[1,1], y[1,1]]   f.dimensions = (time, x, y)

u's time[0,1] isn't native to u at all — it's t's (the modulo
stepping dimension) own [0,1], relabeled onto its parent time by
intervals.promote(lambda d: d.is_SubIterator). t's +1 is only ever
safe because t wraps around a 2-slot circular buffer; it says nothing
about the raw time axis. usave, on the other hand, is genuinely and
exactly defined over time — [0,0] is its real requirement.

The union combines the two into time[0,1]: the harmless promoted offset
inflates the exact save-mode bound. That 1 then flows into
Dimension._arg_values, which subtracts it from usave's own (already
correct) time_M candidate:

loc_maxv = args.get(self.max_name, defaults[self.max_name])  # 4, from usave's own _arg_defaults (save_size - 1)
loc_maxv -= max(interval.upper, 0)                            # 4 - 1 = 3

This fix vs. #2933

#2933 computes the (buggy) union first, then loops back over every
(function, interval-group) pair a second time, intersects each
function's own interval against the already-polluted global one, and
mutates the global interval's upper bound via a new Interval.set_upper/
IntervalGroup.set_upper method introduced for the purpose.

This PR instead prevents the bad union from happening in the first place,
with no new API surface on Interval/IntervalGroup at all — the whole
fix is local to Cluster.dspace:

1. Partition parts into "native" (dimension is genuinely one of the
   function's own .dimensions) vs. the rest.
2. Union only the native parts — this is usave's true time[0,0],
   uncontaminated by u's promoted [0,1].
3. For every Dimension that has a native definer, rebuild its Interval
   keeping the normal, full-union lower but overriding upper with the
   native-only union's upper.

natives = {f: IntervalGroup([i for i in v if i.dim in f.dimensions],
                            relations=v.relations, mode=v.mode)
           for f, v in parts.items()}
natives = {f: v for f, v in natives.items() if v}
if natives:
    native_intervals = IntervalGroup.generate('union', *natives.values())
    rebuilt = [
        Interval(i.dim, i.lower, native_intervals[i.dim].upper, i.stamp)
        if i.dim in native_intervals else i
        for i in intervals
    ]
    intervals = IntervalGroup(rebuilt, relations=intervals.relations,
                              mode=intervals.mode)

Why upper-only isn't an arbitrary special case

I tried it symmetric first (native-only for both bounds), which broke the
existing test_operator.py::TestInternals::test_indirection: dropping a
non-native [0,0] lower contribution shifted dspace[time].lower from 0
to 1.

Dimension._arg_values is already asymmetric between the two bounds:

loc_minv -= min(interval.lower, 0)   # only ever tightens for a *negative* lower
loc_maxv -= max(interval.upper, 0)   # always applies a positive upper

A promoted SubIterator's lower offset is consumed correctly regardless of
promotion (the min(x, 0) clamp makes a non-negative lower a no-op
either way); only its upper offset can wrongly shrink an unrelated
native Function's bound. Restricting the fix to upper matches an
asymmetry that already provably exists in the consumer, rather than being
invented for this fix.

Testing

• Issue default time_M value in presence of saved TimeFunction is wrong #2235 reproducer: fixed (time_M == 4, apply() gives the
  expected data).
• New regression test: test_dimension.py::TestBufferedDimension::test_default_timeM_with_saved.
• tests/test_dimension.py + tests/test_operator.py + tests/test_ir.py:
  429 passed, 0 failed (includes test_indirection, which a symmetric
  version of the fix broke).
• tests/test_adjoint.py: 74 passed, 0 failed.
• Broader sweep (test_dse, test_buffering, test_checkpointing,
  test_lower_clusters, test_lower_exprs, test_dle, test_visitors,
  test_subdomains, test_interpolation): 1004 passed, 2 pre-existing
  xfailed, 0 failed.
• 1507 tests passed total, 0 failures.

[FabioLuporini]

This line is where this comment starts becoming hard to understand
IMO, this whole comment could be replaced by a strightforward example, such as

# Handle inconsistencies due to nonsaved- and saved-TimeFunctions, if any
# E.g., given `f(t, x, y)`  and `fsave(time, x, y)`, ...
# ...

should be a 3~4 lines long comment with an example to stand the change of getting understood. Otherwise, a simple :

# Handle inconsistencies due to nonsaved- and saved-TimeFunctions, if any; see issue #2235

[FabioLuporini]

likely just

key = lambda ...
natives = {f: intervals.project(key) for f, intervals in parts.items()}

[FabioLuporini]

OK so I got to this point and I have a very simple question; if there are both nonsaved- and saved-TimeFunction, instead of this pretty complicated logic, why don't you simply filter the nonsaved TimeFunctions off parts before constructing the intervals? something along the lines of

from devito.tools import split
....
....
# <comment as above>
functions = {f for f in parts if f.is_TimeFunction}
nonsaved, saved = split(functions, lambda f: f.save is None)
if saved:
    parts = {f: v for f, v in parts.items() if f not in nonsaved}