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?

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?

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?

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

EdCaunt · 2025-07-07T08:52:04Z

devito/timestepping/superstep.py

+
+def superstep_generator_iterative(field, stencil, k, tn=0):
+    ''' Generate superstep iteratively:
+    A^j+1 = A·A^j


Some brackets here would improve readability and reduce ambiguity

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)

EdCaunt · 2025-07-07T09:00:37Z

devito/timestepping/superstep.py

+    return u, u_prev, Eq(u.forward, stencil_next), Eq(u_prev.forward, current)
+
+
+def superstep_generator(field, stencil, k, tn=0):


Why is there both superstep_generator and superstep_generator_iterative? There is also some code overlap between these two functions which should be consolidated. Perhaps a SuperStepGenerator class with suitable methods?

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

devito/timestepping/superstep.py

+    ss_stencil = stencil.subs({field: u, field.backward: u_prev}, postprocess=False)
+    ss_stencil = ss_stencil.expand().expand(add=True, nest=True)
+
+    # Binary decomposition algorithm


Nitpick: this could use more comments

EdCaunt · 2025-07-07T09:05:21Z

devito/timestepping/superstep.py

+    ''' Transfer state from a previous TimeFunction to a 2 field superstep
+    Used after injecting source using standard timestepping.
+    '''
+    idx = tn % 3 if old.save is None else -1


This should check if there is a buffer, then either check the buffer size or take the last index, rather than hardcoding 3

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?

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

mloubout · 2025-07-07T13:08:44Z

devito/timestepping/superstep.py

+
+
+def superstep_generator_iterative(field, stencil, k, tn=0):
+    ''' Generate superstep iteratively:


Nitpicking: For consistency with codebase, use """ and start doctrine on next line"

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)

mloubout · 2025-07-07T15:16:05Z

devito/timestepping/superstep.py

+
+    # 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)


why double expand?

mloubout · 2025-07-07T15:16:46Z

devito/timestepping/superstep.py

+                {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)


mloubout · 2025-07-07T15:17:01Z

devito/timestepping/superstep.py

+
+    if k >= 2:
+        for _ in range(k - 2):
+            current = current.subs(


split in two lines for readability?

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

mloubout · 2025-07-07T15:23:11Z

examples/timestepping/acoustic_superstep_2d.py

+    # Stencil and operator
+    idx = tn1 % 3
+    if step == 1:
+        # Non-superstep case


I don't get that case, it just copies into a new Function to do the same thing?

mloubout · 2025-07-07T15:23:42Z

examples/timestepping/acoustic_superstep_2d.py

+    else:
+        new_u, new_u_p, *stencil = superstep_generator(u, stencil.rhs, step)
+
+        new_u.data[0, :] = u.data[idx - 1]


why isn't superstep_generator doing that?

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:

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

devito/timestepping/superstep.py

+    b_tmp = Function(name="b_tmp", grid=grid, space_order=2*k)
+
+    if k >= 2:
+        for _ in range(k - 2):


[nitpick] The loop variable is unused. Consider using a more descriptive name like 'step' or 'iteration' to improve code readability, or use 'for step in range(k - 2):' if the index might be useful for debugging.

Suggested change

for _ in range(k - 2):

for step in range(k - 2):

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

devito/timestepping/superstep.py

+    # 3 should be replaced with `old.time_order + 1` although this needs some thought
+    idx = tn % 3 if old.save is None else -1


This TODO comment indicates incomplete implementation. The hardcoded value '3' should be replaced with a proper calculation based on time_order, or this limitation should be documented if intentional.

Suggested change

# 3 should be replaced with `old.time_order + 1` although this needs some thought

idx = tn % 3 if old.save is None else -1

# Dynamically calculate the modulo base using `old.time_order + 1`

idx = tn % (old.time_order + 1) if old.save is None else -1

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

examples/timestepping/acoustic_superstep_2d.py

+
+def marmousi(grid, step=0):
+    # Grab dataset from pwd or download
+    url = 'https://github.com/devitocodes/data/raw/refs/heads/master/Simple2D/vp_marmousi_bi'  # noqa: E501


[nitpick] The long URL should be broken into multiple lines or stored as a constant to improve readability instead of using noqa comment.

Suggested change

url = 'https://github.com/devitocodes/data/raw/refs/heads/master/Simple2D/vp_marmousi_bi' # noqa: E501

url = MARMOUSI_URL

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]]))

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

Codecov Report

✅ All modified and coverable lines are covered by tests.
⚠️ Please upload report for BASE (main@dd1e6c7). Learn more about missing BASE report.
⚠️ Report is 19 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2658   +/-   ##
=======================================
  Coverage        ?   83.76%           
=======================================
  Files           ?      248           
  Lines           ?    49625           
  Branches        ?     4366           
=======================================
  Hits            ?    41569           
  Misses          ?     7274           
  Partials        ?      782

Flag	Coverage Δ
pytest-gpu-nvc-nvidiaX	`72.07% <ø> (?)`

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?

FabioLuporini · 2025-08-11T08:39:42Z

devito/timestepping/superstep.py

+
+
+def superstep_generator_iterative(field, stencil, k, tn=0):
+    ''' Generate superstep iteratively:


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

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

examples/timestepping/acoustic_superstep_2d.py

+
+import matplotlib.pyplot as plt
+import numpy as np
+from devito import (


blankline between third-party imports and devito imports

FabioLuporini · 2025-08-11T08:45:03Z

examples/timestepping/acoustic_superstep_2d.py

+    pde = (1/velocity**2)*u.dt2 - u.laplace
+    stencil = Eq(u.forward, solve(pde, u.forward))
+
+    tn1 = int(np.ceil((t1 - t0)/model.critical_dt))


nitpicking, for homogeneity: why not just tn? rather, "nt", since it's "number of tiemsteps"?

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 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 added 8 commits August 4, 2025 14:45

dsl: Add superstep generators to the dsl

7051473

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

9bec07b

examples: Add a simple 2D superstep example

6f7dffe

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

4a14987

misc: Init

6e7377f

dsl: Add superstep_solution_transfer function + tweaks

21724e6

examples: Add the first draft of the superstepping notebook

e29d36a

lint: Fix my linter

fb72ca2

JDBetteridge and others added 3 commits August 4, 2025 14:45

misc: Noodling

d8fabb4

WIP

e69f9ba

misc: conflict?

b1b0846

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

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


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



		def superstep_generator_iterative(field, stencil, k, tn=0):
		''' Generate superstep iteratively:



		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

		# 3 should be replaced with `old.time_order + 1` although this needs some thought
		idx = tn % 3 if old.save is None else -1

	url = 'https://github.com/devitocodes/data/raw/refs/heads/master/Simple2D/vp_marmousi_bi' # noqa: E501
	url = MARMOUSI_URL

	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

Are you sure you want to change the base?

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!

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!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull Request Overview

Reviewed Changes

Uh oh!

Copilot AI Jul 24, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Jul 24, 2025

Choose a reason for hiding this comment

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