Skip to content

api: Fix handling of staggering in injection/interpolation#2936

Open
mloubout wants to merge 3 commits into
mainfrom
inject-interp
Open

api: Fix handling of staggering in injection/interpolation#2936
mloubout wants to merge 3 commits into
mainfrom
inject-interp

Conversation

@mloubout
Copy link
Copy Markdown
Contributor

No description provided.

@mloubout mloubout added the API api (symbolics, types, ...) label May 20, 2026
@codecov
Copy link
Copy Markdown

codecov Bot commented May 20, 2026

Codecov Report

❌ Patch coverage is 99.18831% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 83.39%. Comparing base (68e2f74) to head (551d050).

Files with missing lines Patch % Lines
devito/finite_differences/derivative.py 0.00% 1 Missing and 1 partial ⚠️
tests/test_interpolation.py 99.61% 1 Missing and 1 partial ⚠️
devito/types/sparse.py 96.55% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #2936      +/-   ##
==========================================
+ Coverage   83.33%   83.39%   +0.06%     
==========================================
  Files         248      248              
  Lines       51799    51881      +82     
  Branches     4467     4479      +12     
==========================================
+ Hits        43165    43267     +102     
+ Misses       7882     7861      -21     
- Partials      752      753       +1     
Flag Coverage Δ
pytest-gpu-aomp-amdgpuX 68.70% <88.04%> (+0.03%) ⬆️
pytest-gpu-gcc- 78.09% <99.18%> (+0.06%) ⬆️
pytest-gpu-icx- 78.00% <99.18%> (+0.05%) ⬆️
pytest-gpu-nvc-nvidiaX 69.23% <88.04%> (+0.03%) ⬆️

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
Copy link
Copy Markdown

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@mloubout mloubout force-pushed the inject-interp branch 3 times, most recently from bc601dc to 74db52a Compare May 21, 2026 03:08
Comment thread devito/operations/interpolators.py Outdated

def _field_shifts(self, field):
"""
Per-grid-Dimension half-cell shift induced by ``field``'s staggering
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ultra-nitpick, fixable in a separate PR -- single `, not two, for homogeneity

Copy link
Copy Markdown
Contributor

@EdCaunt EdCaunt May 21, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We seem to be wildly inconsistent on this

Comment thread devito/operations/interpolators.py Outdated
try:
staggered = field.staggered
except AttributeError:
return None
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

when can u end up in the AttributeError? add a comment?

Comment thread devito/operations/interpolators.py Outdated
try:
staggered = field.staggered
except AttributeError:
return None
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

return empty tuple

Comment thread devito/operations/interpolators.py Outdated
except AttributeError:
return None
if not staggered or all(s == 0 for s in staggered):
return None
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

return empty tuple

Comment thread devito/operations/interpolators.py Outdated
for _, g in groupby(pairs, lambda f: f[0].staggered):
g = list(g)
g_fields = [f for f, _ in g]
g_exprs = [e for _, e in g]
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you can do

g_fields, g_exprs = zip(*g)

instead of these two lines

Comment thread devito/types/sparse.py Outdated
def _grid_map(self):
return {}

@property
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cached_property

Comment thread devito/types/sparse.py

@property
def origin(self):
return DimensionTuple(*[0]*len(self.dimensions),
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ultra-nitpick: self.ndim == len(self.dimensions)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No this is a SparseFunction, self.ndim is the number of grid dimensions, this is the number of sdimension of the SparseFunction

Comment thread devito/types/sparse.py
def _pos_symbols(self):
return [Symbol(name=f'pos{d}', dtype=np.int32)
for d in self.grid.dimensions]
@memoized_meth
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the fact that this had to be turned into a memoized_meth makes me think that it should better be placed in the interpolator classes instead.
Accepting shifts here doesn't make much sense imho.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I probably agree with this

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the symbol for the sparse position. And it's not used only for interpolation, it's used for plain guards such as here that are used for guarding generic sparse expressions based on the position

Comment thread devito/types/sparse.py
@cached_property
def _position_map(self):
@memoized_meth
def _position_map(self, shifts=None):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same story, imho this method doesn't belong here anymore

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It does, this is a SparseFunction information which is it's position relative to an origin. The interpolator uses it, it doesn't define it.

"assert np.isclose(norm(rec), 31.23, atol=0, rtol=1e-3)\n",
"assert np.isclose(norm(rec2), 3.5482, atol=0, rtol=1e-3)\n",
"assert np.isclose(norm(rec3), 4.7007, atol=0, rtol=1e-3)"
"assert np.isclose(norm(rec), 29.538, atol=0, rtol=1e-3)\n",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that's quite the change

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This one inject the source into v_x which changes things quite a bit over long propagation (2 seconds here) since the source position is not properly nterpolated into the staggered location of v_x

Comment thread devito/types/sparse.py

@property
def origin(self):
return DimensionTuple(*[0]*len(self.dimensions),
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ultra nitpick: DimensionTuple(*[0 for _ in self.dimensions]) would be more readable in this context

Comment thread devito/types/sparse.py
def _pos_symbols(self):
return [Symbol(name=f'pos{d}', dtype=np.int32)
for d in self.grid.dimensions]
@memoized_meth
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I probably agree with this

Comment thread devito/types/sparse.py

# Position map and temporaries for it
pmap = self._position_map
pmap = self._position_map()
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this never need to know about the shifts then?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a guard for pure sparse operation, there is no staggered origin involved

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How come this file has been so heavily reworked?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just reorganized it was a pain to navigate it, it's 99% the same just lifted in classes

'(x, y)', '(x, z)', '(y, z)',
'(x, y, z)'
])
def test_interpolate_staggered(self, stagg):
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note: These three tests are the only new ones in this file

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

API api (symbolics, types, ...)

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants