api: fix staggered evaluation for add/mul #2782

mloubout · 2025-10-31T14:21:29Z

No description provided.

codecov · 2025-10-31T14:41:58Z

Codecov Report

❌ Patch coverage is 76.02740% with 35 lines in your changes missing coverage. Please review.
✅ Project coverage is 83.01%. Comparing base (d88639b) to head (26490d4).
⚠️ Report is 9 commits behind head on main.

Files with missing lines	Patch %	Lines
devito/arch/archinfo.py	18.75%	13 Missing ⚠️
devito/passes/clusters/implicit.py	0.00%	10 Missing ⚠️
devito/types/dense.py	71.42%	3 Missing and 1 partial ⚠️
devito/finite_differences/differentiable.py	72.72%	2 Missing and 1 partial ⚠️
devito/passes/iet/misc.py	78.57%	3 Missing ⚠️
devito/passes/iet/parpragma.py	0.00%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2782      +/-   ##
==========================================
- Coverage   83.03%   83.01%   -0.03%     
==========================================
  Files         248      248              
  Lines       50898    50957      +59     
  Branches     4485     4489       +4     
==========================================
+ Hits        42265    42303      +38     
- Misses       7865     7883      +18     
- Partials      768      771       +3

Flag	Coverage Δ
pytest-gpu-aomp-amdgpuX	`68.49% <54.12%> (-0.05%)`	⬇️
pytest-gpu-nvc-nvidiaX	`69.02% <54.12%> (-0.04%)`	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

review-notebook-app · 2025-11-01T02:55:53Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

FabioLuporini

Uncontroversial

JDBetteridge

A lot of the norms in the tests and notebooks have changed significantly, how do we know these are now the "right values". And is the test really valuable if we do not know what the true value is?

One could suggest that these are in some sense "regression tests", but then we definitely should not be changing the values!

mloubout · 2025-11-06T03:25:56Z

A lot of the norms in the tests and notebooks have changed significantly

The one norm that changed a lot in test_mpi.py .... was indeed because of a bug introduced, made me find it so thanks, the rest is small changes due to a few extra interpolation but within tolerance

FabioLuporini · 2025-11-06T08:12:24Z

devito/finite_differences/differentiable.py



 class DiffDerivative(IndexDerivative, DifferentiableOp):
    pass


drop pass now

FabioLuporini · 2025-11-06T08:13:03Z

devito/finite_differences/differentiable.py

 class DiffDerivative(IndexDerivative, DifferentiableOp):
    pass

+    def _eval_at(self, func):


isn't it the case of IndexDerivative too?

FabioLuporini · 2025-11-06T08:13:33Z

devito/finite_differences/differentiable.py

        kwargs.pop('is_commutative', None)
        return self.func(*args, **kwargs)

+    def _eval_at(self, func):


nitpicking: Perhaps could be moved into a mixin class from which all *Derivative class inherit

Seconded: EvaluatedMixin perhaps?

FabioLuporini · 2025-11-06T08:14:18Z

devito/types/basic.py

            # Two indices are aligned if they differ by an Integer*spacing.
-            v = (i - j)/d.spacing
+            if not i.has(d):
+                # Maybe a subdimension


nitpicking, SubDimension

FabioLuporini · 2025-11-06T08:14:30Z

devito/types/basic.py

-            v = (i - j)/d.spacing
+            if not i.has(d):
+                # Maybe a subdimension
+                dims = {sd for sd in i.free_symbols if getattr(sd, 'is_Dimension', False)


avoids circular (so local) import using getattr

it's just a hack to avoid a hack in that case?

anyway, u just reminded me that unless i is a Number or a Dimension w/o an offset, it will be of type AffineIndexAccessFunction https://github.com/devitocodes/devito/blob/main/devito/types/dimension.py#L1822

so all you have to here is

try: sds, = i.sds except (AttributeError, ValueError): # Expected exactly one StencilDimension continue

btw should we not raise a NotImplementedError if len(sds) > 1 instead?

Don't think it's better. You have to check all three of i.sds, i.d, o.ofs and see if one is a sbdimensions it doesn't make it cleaner or simpler.

btw should we not raise a NotImplementedError if len(sds) > 1 instead?

No if there is multiple dimensions it's considered that the user decided to specifiy its own indices and to be valid ones. Would lead to a lot of breakign cases to raise and error here

Could you use is_Sub?

Could you use is_Sub?
No because it could be something other than a Dimension

FabioLuporini · 2025-11-06T08:17:55Z

devito/types/basic.py

        for the compiler.
        """
+        mapper = self._grid_map
+        subs = mapper.pop('subs', {})


since _grid_map is (must be) a cached_property, use .get() instead of .pop()

In fact, _grid_map, now that you make me think twice about it, should return a frozendict to be pedantic

I should not be cached no, that would break all of the FD since f and f._subs(x, x+.5) would end up with the same _grid_map when the second needs interpolation

using pop so that can just iterate mapper below

oh that's because we still have that crazy hack that makes every instance of f(x,y,z) reference the same __dict__ as the "origin" Function?

No this is because this is based on indices. Every f.subs is the same f with different indices (and has too or it would get really nasty with data/.....) and this checks if those specific indices are on the grid given the reference index (staggered "dim")

FabioLuporini · 2025-11-06T08:18:33Z

devito/types/basic.py

        # Evaluate. Since we used `self.function` it will be on the grid when
        # evaluate is called again within FD
        retval = retval._evaluate(**kwargs)
+        if subs:


I don't think this if is needed, .subs should take O(1) when subs is empty by returning immediately

FabioLuporini · 2025-11-06T08:19:32Z

devito/types/dense.py

+            except KeyError:
+                pass
+
        if mapper:


nitpicking, but same as before. We can reduce verbosity by avoiding these if since they practically speaking won't improve performance

FabioLuporini · 2025-11-06T08:19:43Z

devito/types/dense.py

+
        if mapper:
            return self.subs(mapper)
        return self


nitpicking, perhaps a blank line

FabioLuporini · 2025-11-06T08:20:47Z

tests/test_derivatives.py

        eqne = eqn.evaluate.rhs
        assert simplify(eqne - (p._subs(y, yp).evaluate * f).dx(x0=xp).evaluate) == 0

+    def test_param_stagg_add(self):


quite a few changes have been made, and it seems weird that only one new test was required?

were they all fixing the same exact (set of) thing(s), which is entirely captured by this test? if so, OK

This is the only "change" that is tested. Everything else isadjustments from removing is_Parameter and using subs in eval_for_fd instead of explicit derivatives

JDBetteridge

Just flagging a couple of other big changes that might need investigation (or given they are just in the notebooks, maybe re-running after your bug fix)

JDBetteridge · 2025-11-06T12:15:08Z

examples/cfd/02_convection_nonlinear.ipynb

   "source": [
    "from devito import norm\n",
-    "assert np.isclose(norm(u), norm(U[0]), rtol=1e-5, atol=0)"
+    "assert np.isclose(norm(u), norm(U[0]), rtol=1e-2, atol=0)"


This is 3 orders of magnitude change in the difference between two values that are notionally the same

The numpy implementation is techinically wrong, as it does staggering but does not interpolate u/v where needed so there is some difference now because forcing the devito equation to be the same would be tricky since "incorrect"

JDBetteridge · 2025-11-06T12:15:39Z

examples/userapi/07_functions_on_subdomains.ipynb

   "outputs": [],
   "source": [
-    "assert np.isclose(np.linalg.norm(rec.data), 4263.511, atol=0, rtol=1e-4)"
+    "assert np.isclose(np.linalg.norm(rec.data), 3640.584, atol=0, rtol=1e-4)"


This jump is also pretty big!

Likely some parameters interpolation that were not being processed due to the SubDimensions. Will double check

Yes, basically comes from parameters (damp, ro, ...) being properly interpolated now.

…rtree at read

devito/finite_differences/differentiable.py

FabioLuporini · 2025-12-02T09:49:32Z

devito/passes/clusters/implicit.py

+    def key(i):
+        try:
+            return i.indices[d]
+        except (KeyError, AttributeError):


why? and what triggers this?

Same as above

FabioLuporini · 2025-12-02T09:50:21Z

devito/passes/clusters/implicit.py

+                edims = set(retrieve_dimensions(tkn, deep=True))
+                return dim._defines & edims and edims.issubset(prefix.dimensions)
+
+            mapper = {k: v for k, v in mapper.items() if key(v)}


what triggered this change?

Current version is too simplistic and misses things like combination of implicit dims. Those changes here are preleminary to the topography multisubdomain when for example taken in intersection of an abox where you end up with multimple implicit dims in the tkn.

FabioLuporini · 2025-12-02T09:52:36Z

devito/passes/iet/misc.py

+    except (AttributeError, ValueError):
+        dtype = np.float32
+        eps = np.finfo(dtype).resolution**2
+    b = printer()._print(Cast('b', dtype=dtype))


can use a blank line before this line

also, instantiating the printer at this point is quite weird. Even more if you think one new instance will be instantiated each time you enter this function. Why don't we create the printer once and for all at the beginning of lowering or wherever (as soon as ...) possible?

That's how the printer is used throughout the codebase. It's instantiated when needed.

devito/passes/iet/misc.py

FabioLuporini · 2025-12-02T09:59:41Z

devito/types/basic.py

+        retval = self.subs({i.subs(subs): self.indices_ref[d]
+                            for d, i in mapper.items()})
+        if 'harmonic' in self._avg_mode:
+            retval = retval.safe_inv(retval, safe='safe' in self._avg_mode)


safe='safe' is a bit clumsy

It's a boolean 'safe' in self._avg_mode not safe='safe'

FabioLuporini · 2025-12-02T10:01:03Z

devito/types/basic.py

-        if self._avg_mode == 'harmonic':
-            from devito.finite_differences.differentiable import SafeInv
-            retval = SafeInv(retval, self.function)
+        if 'harmonic' in self._avg_mode:


nitpicking: the _avg_mode should be an instance of something rather than an iterable of string, so that you may do e.g. self._avg_mode.is_harmonic or self._avg_mode.is_harmonic_safe

FabioLuporini · 2025-12-02T10:05:00Z

devito/types/dense.py

-            return self.subs(mapper)
-        return self
+
+        mapper = {}


aka

m0 = self.indices_ref m1 = func.indices_ref mapper = {m[d]: m1[d] for d in self.dimensions if m0.get(d) is not m1.get(d)}

?

No that would still throw a KeyError when d is in m0 and not m1 as m0.get(d) is not m1.get(d) would be True but m1[d] doesn't exist

FabioLuporini · 2025-12-02T10:05:54Z

devito/types/dense.py

+        retval = super()._evaluate(**kwargs)
+        if not self._time_buffering and not retval.is_Function:
+            # Saved TimeFunction might need streaming, expand interpolations
+            # for easier processing.


(ultra nitpick) no full stop at the end of comments

FabioLuporini · 2025-12-02T10:06:06Z

devito/types/dense.py

+    def _evaluate(self, **kwargs):
+        retval = super()._evaluate(**kwargs)
+        if not self._time_buffering and not retval.is_Function:
+            # Saved TimeFunction might need streaming, expand interpolations


why the retval.is_Function part? what else can it be?

Derivative when there is an interp. But if it's a Function then you get an infinite recursion on .evaluate

EdCaunt · 2025-12-05T15:44:10Z

devito/arch/archinfo.py

+        lib = ctypes.CDLL(libname)
+
+        cudaDevAttrMaxSharedMemoryPerBlockOptin = 97
+        # get current device


Ultra nitpick: Comment doesn't start with a capital

EdCaunt · 2025-12-05T15:44:58Z

devito/arch/archinfo.py

+            return 64 * 1024  # 64 KB default
+        lib = ctypes.CDLL(libname)
+
+        cudaDevAttrMaxSharedMemoryPerBlockOptin = 97


Maybe abbreviate this name and add a comment?

That's the actual CUDA name

EdCaunt · 2025-12-05T15:48:28Z

devito/finite_differences/differentiable.py

        kwargs.pop('is_commutative', None)
        return self.func(*args, **kwargs)

+    def _eval_at(self, func):


Seconded: EvaluatedMixin perhaps?

EdCaunt · 2025-12-05T15:57:38Z

devito/types/basic.py

-        if self._avg_mode not in ['arithmetic', 'harmonic']:
+        if self._avg_mode not in ['arithmetic', 'harmonic', 'safe_harmonic']:
            raise ValueError("Invalid averaging mode_mode %s, accepted values are"
                             " arithmetic or harmonic" % self._avg_mode)


Does this error need updating?

EdCaunt · 2025-12-05T16:05:46Z

devito/types/basic.py

-            v = (i - j)/d.spacing
+            if not i.has(d):
+                # Maybe a subdimension
+                dims = {sd for sd in i.free_symbols if getattr(sd, 'is_Dimension', False)


Could you use is_Sub?

mloubout added the API api (symbolics, types, ...) label Oct 31, 2025

mloubout force-pushed the fd-eval-add branch from f709216 to be8bdbf Compare October 31, 2025 14:25

mloubout force-pushed the fd-eval-add branch 2 times, most recently from 766cc81 to a156d3c Compare October 31, 2025 16:05

mloubout force-pushed the fd-eval-add branch 2 times, most recently from 1a9951a to dc0e8f5 Compare November 3, 2025 13:01

FabioLuporini approved these changes Nov 3, 2025

View reviewed changes

mloubout force-pushed the fd-eval-add branch 7 times, most recently from 42588ae to 6c6d45f Compare November 4, 2025 04:17

JDBetteridge reviewed Nov 5, 2025

View reviewed changes

mloubout force-pushed the fd-eval-add branch from 6c6d45f to 45bc728 Compare November 6, 2025 03:24

FabioLuporini reviewed Nov 6, 2025

View reviewed changes

JDBetteridge requested changes Nov 6, 2025

View reviewed changes

mloubout force-pushed the fd-eval-add branch 5 times, most recently from cfcadc9 to ef1f7d8 Compare November 12, 2025 15:06

mloubout force-pushed the fd-eval-add branch 4 times, most recently from 01d3d1a to cea48f3 Compare November 18, 2025 02:08

mloubout force-pushed the fd-eval-add branch 7 times, most recently from 0823c5f to 7821fc1 Compare November 21, 2025 14:48

mloubout force-pushed the fd-eval-add branch 3 times, most recently from 3cd6f39 to f91130f Compare November 26, 2025 14:06

mloubout added 6 commits December 1, 2025 07:08

api: fix staggered evaluation for add/mul

9077a30

api: remove obsolete and clunky is_parameter

325cca1

examples: remove usage of obsolete paramter kw

71bfb4f

examples: update notebooks

e4b71c5

api: prevent evaluated derivatives to be re-evaluted

7c4f88e

compiler: prevent hosted per-thread arrays are dereferenced within pa…

5b93b20

…rtree at read

mloubout force-pushed the fd-eval-add branch from f91130f to 136f771 Compare December 1, 2025 12:08

FabioLuporini reviewed Dec 2, 2025

View reviewed changes

compiler: patch multi-subdim isapce handling

55f84b2

mloubout force-pushed the fd-eval-add branch 3 times, most recently from b1f1e14 to edd10f2 Compare December 4, 2025 19:44

arch: get shm max size on device

26490d4

mloubout force-pushed the fd-eval-add branch from edd10f2 to 26490d4 Compare December 4, 2025 20:24

EdCaunt reviewed Dec 5, 2025

View reviewed changes

mloubout merged commit a00cca8 into main Dec 8, 2025
39 checks passed

mloubout deleted the fd-eval-add branch December 8, 2025 17:37



		class DiffDerivative(IndexDerivative, DifferentiableOp):
		pass

api: fix staggered evaluation for add/mul #2782

api: fix staggered evaluation for add/mul #2782

Uh oh!

Conversation

mloubout commented Oct 31, 2025

Uh oh!

codecov bot commented Oct 31, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

review-notebook-app bot commented Nov 1, 2025

Uh oh!

FabioLuporini left a comment

Choose a reason for hiding this comment

Uh oh!

JDBetteridge left a comment

Choose a reason for hiding this comment

Uh oh!

mloubout commented Nov 6, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JDBetteridge left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Oct 31, 2025 •

edited

Loading

mloubout Dec 2, 2025 •

edited

Loading