dsl: Superstep #2658

JDBetteridge · 2025-07-03T15:15:31Z

Currently a PR into #2559 as the diff is smaller

Adds superstep functionality as well as example scripts and notebook

review-notebook-app · 2025-07-03T15:15:37Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

review-notebook-app · 2025-07-04T10:14:01Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:01Z
----------------------------------------------------------------

"In Nemeth et al k=30 is chosen."

review-notebook-app · 2025-07-04T10:14:02Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:01Z
----------------------------------------------------------------

Line #3.    src_term = source.inject(field=u.forward, expr=source*velocity*velocity*dt*dt)

expr=source*velocity**2*dt**2

review-notebook-app · 2025-07-04T10:14:02Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:02Z
----------------------------------------------------------------

Typo: arguemnts

"the size of the superstep: k in the paper, superstep_size in the code" <- add a colon

"Written in maths:" <- odd phrasing, reconsider

review-notebook-app · 2025-07-04T10:14:03Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:03Z
----------------------------------------------------------------

Line #1.    # Should the timer context be removed???

Leftover?

JDBetteridge commented on 2025-08-19T12:37:51Z
----------------------------------------------------------------

More a question for the reviewer, I guess the answer is: Yes, remove it

review-notebook-app · 2025-07-04T10:14:04Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:03Z
----------------------------------------------------------------

Line #7.    # Fetch and setup the

Set up the what?

JDBetteridge commented on 2025-08-19T12:39:44Z
----------------------------------------------------------------

yes

EdCaunt commented on 2025-08-20T10:48:20Z
----------------------------------------------------------------

This comment is still missing a noun

review-notebook-app · 2025-07-04T10:14:04Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:04Z
----------------------------------------------------------------

"we would have to redefine much of what we have done above."

What would need redefining?

JDBetteridge commented on 2025-08-19T12:42:21Z
----------------------------------------------------------------

I'll double check, I'm fairly sure something was unstable

JDBetteridge commented on 2025-08-19T12:57:03Z
----------------------------------------------------------------

Apparently nothing, I'll update to use a different dt

EdCaunt commented on 2025-08-20T10:49:30Z
----------------------------------------------------------------

Note: this is still in the notebook

review-notebook-app · 2025-07-04T10:14:05Z

View / edit / reply to this conversation on ReviewNB

EdCaunt commented on 2025-07-04T10:14:05Z
----------------------------------------------------------------

x10^-4

JDBetteridge commented on 2025-08-19T12:43:29Z
----------------------------------------------------------------

Technically 10^4 isn't wrong...

EdCaunt commented on 2025-08-20T10:52:35Z
----------------------------------------------------------------

It is misleading though and not how most people would interpret that statement

devito/timestepping/superstep.py

EdCaunt · 2025-07-07T08:58:25Z

devito/timestepping/superstep.py

+    # New fields, for vector formulation both current and previous timestep are needed
+    name = field.name
+    grid = field.grid
+    u = TimeFunction(name=f'{name}_ss', grid=grid, time_order=2, space_order=2*k)


I presume the space_order=2*k is to ensure the halo is sufficiently wide for the operators? Would taking the space order from the original fields and using padding here be strictly more correct?

This needs to use ._rebuild (or .func) to make sure it's consistent with field

u = field._rebuild(name=f'{name}_ss', space_order=2*k) u_prev = field._rebuild(name=f'{name}_ss', space_order=2*k)

So, something is fundamentally different about doing it the way you suggest @mloubout.

Previously:

With the suggested change:

Seems to be an off by k error

Yeah actually ran into that in something else. _rebuild keeps halo from the original func so the order is right but the halo is too short so it leads to issues unless the new halo is specified as well.

devito/timestepping/superstep.py

EdCaunt · 2025-07-07T09:06:12Z

devito/timestepping/superstep.py

+    '''
+    # Placeholder fields for forming the superstep
+    grid = u.grid
+    a_tmp = Function(name="a_tmp", grid=grid, space_order=2*k)


Are these Functions just used as symbolic objects? If so, why not just use Symbol?

It doesn't work, trying TempFunction although I have no idea if it's the correct type to try

TempFunction also doesn't work, I'm honestly not sure what using a different object gains you here. These temporaries are never used, so should never be allocated any memory. Unless I'm misunderstanding something?

EdCaunt · 2025-07-07T09:08:47Z

examples/timestepping/acoustic_superstep_2d.py

+        data = acoustic_model(model, step=step, snapshots=m)
+        time = np.linspace(model.t1, model.t2, (m - 1)//step + 1)
+        idx = 0
+        for ii, ax in enumerate(ax_row):


I would set up a separate plotting function for the contents of this loop to improve readability

The plotting isn't really the bulk of this, its mainly setting different axis properties for different iteration counts. I can try and refactor it further

devito/timestepping/superstep.py

mloubout · 2025-07-07T13:10:38Z

devito/timestepping/superstep.py

+    # New fields, for vector formulation both current and previous timestep are needed
+    name = field.name
+    grid = field.grid
+    u = TimeFunction(name=f'{name}_ss', grid=grid, time_order=2, space_order=2*k)


This needs to use ._rebuild (or .func) to make sure it's consistent with field

u = field._rebuild(name=f'{name}_ss', space_order=2*k) u_prev = field._rebuild(name=f'{name}_ss', space_order=2*k)

devito/timestepping/superstep.py

mloubout · 2025-07-07T15:21:13Z

examples/timestepping/acoustic_superstep_2d.py

+
+
+def marmousi(grid, step=0):
+    # Grab dataset from pwd or download


This should just throw a warning/erro people should git clone the data repo

And it's a preset in the Model as well so this should be reusing the examples

I might try making a more generic function in examples, but this really needs to go in a separate PR as it's unrelated to this work.

examples/timestepping/acoustic_superstep_2d.py

mloubout · 2025-07-07T15:24:03Z

examples/timestepping/acoustic_superstep_2d.py

+    return u_save.data
+
+
+layered_model = model(


USe existing preset

mloubout · 2025-07-07T15:26:24Z

examples/timestepping/superstep_2d.py

@@ -0,0 +1,114 @@
+''' Script that demonstrates the functionality of the superstep in 2D


That's a lot of duplicate from the 1D case, can probably be merged into one file with 1d/2d as input

Copilot

Pull Request Overview

This PR introduces superstep functionality to Devito, enabling multi-step time advancement with larger effective time steps while maintaining stability. The implementation includes both iterative and binary decomposition algorithms for generating superstep operators.

Adds core superstep generation functionality with two algorithmic approaches
Provides comprehensive examples demonstrating 1D and 2D wave propagation scenarios
Updates test configuration to support higher-order spatial discretizations

Reviewed Changes

Copilot reviewed 5 out of 7 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
devito/timestepping/superstep.py	Core superstep implementation with generator functions and solution transfer utilities
examples/timestepping/superstep_1d.py	1D wave equation example demonstrating superstep functionality
examples/timestepping/superstep_2d.py	2D ripple simulation example with visualization
examples/timestepping/acoustic_superstep_2d.py	Advanced 2D acoustic modeling with source injection and velocity models
tests/test_pickle.py	Updates space_order from 2 to 4 for derivative testing

Comments suppressed due to low confidence (2)

examples/timestepping/superstep_1d.py:49

[nitpick] Variable name 'newu' is inconsistent with the naming convention used elsewhere in the codebase. Consider using 'new_u' to match the pattern used in line 55 and other files.

        newu = u

examples/timestepping/acoustic_superstep_2d.py:28

Class name 'model' should follow Python naming conventions and be capitalized as 'Model'.

class model:

devito/timestepping/superstep.py

examples/timestepping/acoustic_superstep_2d.py

Copilot · 2025-07-24T08:54:21Z

examples/timestepping/acoustic_superstep_2d.py

+                    )))
+                    ax.set_xlim(model.origin[0], model.origin[0] + model.extent[0])
+                    yticks = ax.get_yticks()
+                    ax.set_yticks(np.concat(((model.origin[1],), yticks[2::2])))


[nitpick] Using np.concat with nested tuples is unnecessarily complex. Consider using np.concatenate([model.origin[1:2], yticks[2::2]]) or np.hstack for better readability.

Suggested change

ax.set_yticks(np.concat(((model.origin[1],), yticks[2::2])))

ax.set_yticks(np.concatenate([model.origin[1:2], yticks[2::2]]))

np.concat and np.concatenate are the same function! 🙄

codecov · 2025-08-05T13:31:13Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 82.94%. Comparing base (b90bb2b) to head (a43feb7).
⚠️ Report is 29 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2658   +/-   ##
=======================================
  Coverage   82.94%   82.94%           
=======================================
  Files         248      248           
  Lines       50022    50022           
  Branches     4402     4402           
=======================================
  Hits        41492    41492           
  Misses       7776     7776           
  Partials      754      754

Flag	Coverage Δ
pytest-gpu-aomp-amdgpuX	`68.81% <ø> (+0.01%)`	⬆️
pytest-gpu-nvc-nvidiaX	`69.34% <ø> (ø)`

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.

FabioLuporini

Somewhere, add a docstring referencing Tamas 's paper? or have I missed it?

Perhaps add it to a short README to examples/timestepping, which would also help navigating the directory?

devito/timestepping/superstep.py

FabioLuporini · 2025-08-11T08:44:01Z

examples/timestepping/acoustic_superstep_2d.py

+    solve,
+)
+from devito.timestepping.superstep import superstep_generator
+from scipy.interpolate import interpn


imports from third-parties should go in between those from the stdlin and the devito ones

see here https://github.com/devitocodes/devito/blob/main/CONTRIBUTING.md#coding-guidelines

Cheers.

Any reason that this file isn't linked in either the readme or the wiki?

Can I suggest moving to ruff again, or if we don't like ruff, just using isort? Lots of the existing "Devito style" nitpicks could quite easily be caught by the linter.

examples/timestepping/acoustic_superstep_2d.py

FabioLuporini · 2025-08-11T08:47:36Z

devito/timestepping/superstep.py

+    new.data[0, :] = old.data[idx - 1]
+    new.data[1, :] = old.data[idx]
+    new_p.data[0, :] = old.data[idx - 2]
+    new_p.data[1, :] = old.data[idx - 1]


yes, this is very suspicious...
However, you may as well throw a warning for now

JDBetteridge · 2025-08-18T16:19:29Z

devito/timestepping/superstep.py

+def superstep_generator_iterative(field, stencil, k, tn=0):
+    ''' Generate superstep iteratively:
+    Aʲ⁺¹ = A·Aʲ
+    '''
+    # New fields, for vector formulation both current and previous timestep are needed
+    name = field.name
+    grid = field.grid
+    u = TimeFunction(name=f'{name}_ss', grid=grid, time_order=2, space_order=2*k)
+    u_prev = TimeFunction(name=f'{name}_ss_p', grid=grid, time_order=2, space_order=2*k)
+
+    superstep_solution_transfer(field, u, u_prev, tn)
+
+    # Substitute new fields into stencil
+    ss_stencil = stencil.subs({field: u, field.backward: u_prev}, postprocess=False)
+    ss_stencil = ss_stencil.expand().expand(add=True, nest=True)
+    current = ss_stencil
+
+    # Placeholder fields for forming the superstep
+    a_tmp = Function(name="a_tmp", grid=grid, space_order=2*k)
+    b_tmp = Function(name="b_tmp", grid=grid, space_order=2*k)
+
+    if k >= 2:
+        for _ in range(k - 2):
+            current = current.subs(
+                {u: a_tmp, u_prev: b_tmp}, postprocess=False).subs(
+                {a_tmp: ss_stencil, b_tmp: u}, postprocess=False
+            )
+            current = current.expand().expand(add=True, nest=True)
+    else:
+        current = u
+
+    stencil_next = current.subs(
+        {u: a_tmp, u_prev: b_tmp}, postprocess=False).subs(
+        {a_tmp: ss_stencil, b_tmp: u}, postprocess=False
+    )
+    stencil_next = stencil_next.expand().expand(add=True, nest=True)
+    return u, u_prev, Eq(u.forward, stencil_next), Eq(u_prev.forward, current)
+
+


Suggested change

def superstep_generator_iterative(field, stencil, k, tn=0):

''' Generate superstep iteratively:

Aʲ⁺¹ = A·Aʲ

'''

# New fields, for vector formulation both current and previous timestep are needed

name = field.name

grid = field.grid

u = TimeFunction(name=f'{name}_ss', grid=grid, time_order=2, space_order=2*k)

u_prev = TimeFunction(name=f'{name}_ss_p', grid=grid, time_order=2, space_order=2*k)

superstep_solution_transfer(field, u, u_prev, tn)

# Substitute new fields into stencil

ss_stencil = stencil.subs({field: u, field.backward: u_prev}, postprocess=False)

ss_stencil = ss_stencil.expand().expand(add=True, nest=True)

current = ss_stencil

# Placeholder fields for forming the superstep

a_tmp = Function(name="a_tmp", grid=grid, space_order=2*k)

b_tmp = Function(name="b_tmp", grid=grid, space_order=2*k)

if k >= 2:

for _ in range(k - 2):

current = current.subs(

{u: a_tmp, u_prev: b_tmp}, postprocess=False).subs(

{a_tmp: ss_stencil, b_tmp: u}, postprocess=False

)

current = current.expand().expand(add=True, nest=True)

else:

current = u

stencil_next = current.subs(

{u: a_tmp, u_prev: b_tmp}, postprocess=False).subs(

{a_tmp: ss_stencil, b_tmp: u}, postprocess=False

)

stencil_next = stencil_next.expand().expand(add=True, nest=True)

return u, u_prev, Eq(u.forward, stencil_next), Eq(u_prev.forward, current)

Drop this before merge. I don't think there's much utility to having both methods present

JDBetteridge · 2025-08-18T17:31:29Z

examples/timestepping/ic_superstep.py

I don't think anyone should write code like this and no longer think this is a good example. I may just have to remove it if I can't refactor it into something presentable

JDBetteridge · 2025-08-19T12:37:52Z

More a question for the reviewer, I guess the answer is: Yes, remove it

View entire conversation on ReviewNB

JDBetteridge · 2025-08-19T12:39:45Z

yes

View entire conversation on ReviewNB

JDBetteridge · 2025-08-19T12:42:22Z

I'll double check, I'm fairly sure something was unstable

View entire conversation on ReviewNB

… code

JDBetteridge added documentation API api (symbolics, types, ...) examples examples labels Jul 3, 2025

JDBetteridge requested review from FabioLuporini, ggorman, mloubout and EdCaunt July 3, 2025 15:52

EdCaunt reviewed Jul 7, 2025

View reviewed changes

mloubout reviewed Jul 7, 2025

View reviewed changes

ggorman requested a review from Copilot July 24, 2025 08:53

Copilot AI reviewed Jul 24, 2025

View reviewed changes

JDBetteridge force-pushed the JDBetteridge/expand_derivative branch 2 times, most recently from 4fb9198 to 1268740 Compare July 31, 2025 16:14

Base automatically changed from JDBetteridge/expand_derivative to main July 31, 2025 18:17

JDBetteridge force-pushed the JDBetteridge/superstep branch from 438a601 to b1b0846 Compare August 5, 2025 13:15

FabioLuporini reviewed Aug 11, 2025

View reviewed changes

JDBetteridge commented Aug 18, 2025

View reviewed changes

JDBetteridge and others added 27 commits September 10, 2025 21:34

examples: Add a 1D wave on a string superstep example

90bfd21

examples: Add a simple 2D superstep example

9b6a7de

examples: Add a 2D acoustic superstep example for layered and Marmousi

6149d6f

misc: Init

0170cd4

dsl: Add superstep_solution_transfer function + tweaks

c159970

examples: Add the first draft of the superstepping notebook

80df5de

lint: Fix my linter

9ff9c43

misc: Noodling

b799192

WIP

098f202

misc: conflict?

3a68616

misc: PR comments part 1

7cb4485

misc: Replace more triple quotes

5e3b879

misc: Revert suggestion until fix found

b2fd861

examples: Refactor 1D and 2D code into one file

16bd5d8

examples: Update superstepping script to use the examples.model

e005e6d

misc: Typo

74ea6ea

examples: Review comments on notebook

ce39a89

misc: Add more verbose comment to non-superstepping solution transfer…

b7d4da6

… code

examples: Tidy up 1D/2D example

992abcc

ci: Add superstepping to the CI

11a1d66

misc: Raise exception for time_order not equal 2

a2f478c

misc: Remove dead code

f0b7d99

misc: Add a paper reference to the superstepping module

83d8783

examples: Tidy superstepping notebook

e8496a9

misc: Rename example filename

75eeaf5

ci: Grab correct source from data repo

392e6da

dsl: Remove iterative superstep generator

a43feb7

JDBetteridge force-pushed the JDBetteridge/superstep branch from 28865f3 to a43feb7 Compare September 10, 2025 20:35

mloubout merged commit 0d90d06 into main Sep 11, 2025
36 checks passed

mloubout deleted the JDBetteridge/superstep branch September 11, 2025 12:13



		def marmousi(grid, step=0):
		# Grab dataset from pwd or download

		@@ -0,0 +1,114 @@
		''' Script that demonstrates the functionality of the superstep in 2D

	ax.set_yticks(np.concat(((model.origin[1],), yticks[2::2])))
	ax.set_yticks(np.concatenate([model.origin[1:2], yticks[2::2]]))

dsl: Superstep #2658

dsl: Superstep #2658

Uh oh!

Conversation

JDBetteridge commented Jul 3, 2025

Uh oh!

review-notebook-app bot commented Jul 3, 2025

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app bot commented Jul 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

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!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

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!

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!

Uh oh!

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!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull Request Overview

Reviewed Changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Copilot AI Jul 24, 2025

Choose a reason for hiding this comment

Uh oh!

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

review-notebook-app bot commented Jul 4, 2025 •

edited

Loading

codecov bot commented Aug 5, 2025 •

edited

Loading