26 · Domain Decomposition on Overthrust 3-D¶

A real 3-D benchmark split across GPUs. ModelParallel slices one model into tiles — one per GPU — and exchanges a thin halo every time step, so a single shot is solved cooperatively rather than replicated.

The point of this notebook is not "DD runs" but DD gives the same answer: the record and the gradient are checked against the single-domain result, and the gradient is plotted across the cut lines where a halo bug would show up.

DD needs one process per tile. A notebook is a single process, so the multi-tile part is launched as a subprocess with torch.distributed.run (that is what torchrun is). The world=1 path below runs inline.

In [1]:

import numpy as np, torch, sweep
from sweep.datasets import load
from sweep.parallel import pad_to_mesh

print("GPUs :", torch.cuda.device_count())

# Overthrust 3-D, decimated 2x -> 94 x 401 x 401 @ 50 m (fits one V100 for the
# single-domain reference; DD then splits it py x px).
d  = load("overthrust", variant="3d-acoustic", downsample=2)
vp = np.ascontiguousarray(d["vp"]).astype("float32")
dh = float(d["dh"][0])

# DD v1 slices into UNIFORM tiles, so each split axis must divide evenly by the
# tile count. Overthrust is 401 wide, which 2 and 4 do not divide. Pad -- do
# NOT crop: cropping would silently solve a different, smaller problem.
#
# `pad_to_mesh` pads the HIGH side only (source/receiver indices stay valid)
# and replicates the edge, so the added columns carry no impedance contrast
# and generate no artificial reflection. It is also differentiable: in an
# inversion you keep the UNPADDED model as the optimisation variable and pad
# inside the closure -- gradients on the invented cells flow back onto the
# edge (see test/dd_pad_grad_check.py). Here we pad once for the largest mesh
# factor we intend to use.
TILE = 4
vp_raw = vp
vp = pad_to_mesh(vp, py=TILE, px=TILE)
pad_y, pad_x = vp.shape[1] - vp_raw.shape[1], vp.shape[2] - vp_raw.shape[2]
print("vp", vp_raw.shape, "->", vp.shape,
      f"(edge-padded +{pad_y} in y, +{pad_x} in x for {TILE}-way tiles)")
print(f"   {vp.min():.0f}-{vp.max():.0f} m/s   dh={dh:.0f} m")

import numpy as np, torch, sweep from sweep.datasets import load from sweep.parallel import pad_to_mesh print("GPUs :", torch.cuda.device_count()) # Overthrust 3-D, decimated 2x -> 94 x 401 x 401 @ 50 m (fits one V100 for the # single-domain reference; DD then splits it py x px). d = load("overthrust", variant="3d-acoustic", downsample=2) vp = np.ascontiguousarray(d["vp"]).astype("float32") dh = float(d["dh"][0]) # DD v1 slices into UNIFORM tiles, so each split axis must divide evenly by the # tile count. Overthrust is 401 wide, which 2 and 4 do not divide. Pad -- do # NOT crop: cropping would silently solve a different, smaller problem. # # `pad_to_mesh` pads the HIGH side only (source/receiver indices stay valid) # and replicates the edge, so the added columns carry no impedance contrast # and generate no artificial reflection. It is also differentiable: in an # inversion you keep the UNPADDED model as the optimisation variable and pad # inside the closure -- gradients on the invented cells flow back onto the # edge (see test/dd_pad_grad_check.py). Here we pad once for the largest mesh # factor we intend to use. TILE = 4 vp_raw = vp vp = pad_to_mesh(vp, py=TILE, px=TILE) pad_y, pad_x = vp.shape[1] - vp_raw.shape[1], vp.shape[2] - vp_raw.shape[2] print("vp", vp_raw.shape, "->", vp.shape, f"(edge-padded +{pad_y} in y, +{pad_x} in x for {TILE}-way tiles)") print(f" {vp.min():.0f}-{vp.max():.0f} m/s dh={dh:.0f} m")

GPUs : 4
[sweep.datasets] overthrust:3d-acoustic
  cite:    Aminzadeh, Brac, Kunz (1997), SEG/EAGE 3-D Modeling Series No. 1
  license: CC-BY-4.0

vp (94, 401, 401) -> (94, 404, 404) (edge-padded +3 in y, +3 in x for 4-way tiles)
   2179-6000 m/s   dh=50 m

In [2]:

import matplotlib.pyplot as plt

nz, ny, nx = vp.shape
fig, ax = plt.subplots(1, 3, figsize=(15, 3.4), constrained_layout=True)
for a, (img, ttl, xl, yl) in zip(ax, (
        (vp[nz//3], f"depth slice z={nz//3}", "x", "y"),
        (vp[:, ny//2, :], f"inline y={ny//2}", "x", "z"),
        (vp[:, :, nx//2], f"crossline x={nx//2}", "y", "z"))):
    im = a.imshow(img, cmap="jet", aspect="auto", origin="upper")
    a.set_title(ttl, fontsize=10); a.set_xlabel(xl); a.set_ylabel(yl)
    plt.colorbar(im, ax=a, shrink=0.85, label="vp (m/s)")
plt.show()

import matplotlib.pyplot as plt nz, ny, nx = vp.shape fig, ax = plt.subplots(1, 3, figsize=(15, 3.4), constrained_layout=True) for a, (img, ttl, xl, yl) in zip(ax, ( (vp[nz//3], f"depth slice z={nz//3}", "x", "y"), (vp[:, ny//2, :], f"inline y={ny//2}", "x", "z"), (vp[:, :, nx//2], f"crossline x={nx//2}", "y", "z"))): im = a.imshow(img, cmap="jet", aspect="auto", origin="upper") a.set_title(ttl, fontsize=10); a.set_xlabel(xl); a.set_ylabel(yl) plt.colorbar(im, ax=a, shrink=0.85, label="vp (m/s)") plt.show()

Single-domain reference¶

One GPU, the whole model. This is the answer DD has to reproduce.

In [3]:

from sweep.equations import Acoustic3D
from sweep.propagator.torch import PropTorch
from sweep.propagator.options import BoundarySaving

# 5 Hz Ricker + order-8: fmax (~12.5 Hz) stays above 3.5 ppw at vmin -> no visible grid dispersion
dev, so, abcn, nt, DT = torch.device("cuda:0"), 8, 20, 1500, 0.002

def ricker(nt, dt, fm=5.0, delay=0.25):
    t = np.arange(nt, dtype=np.float32) * dt - delay
    a = np.pi * fm * t
    return ((1.0 - 2.0 * a**2) * np.exp(-(a**2))).astype(np.float32)

wav = torch.as_tensor(ricker(nt, DT), device=dev)
src = np.array([[[nx//2, ny//2, 4]]], dtype=np.int64)              # (x, y, z)
rec = np.array([[[ix, ny//2, 4] for ix in range(6, nx-6, 2)]], dtype=np.int64)

def build(shape):
    return PropTorch(Acoustic3D(spatial_order=so, device=dev, backend="torch"),
                     backend="torch", impl="c", shape=shape, dev=dev, dh=dh, dt=DT,
                     source_type=["h1"], receiver_type=["h1"], abcn=abcn,
                     free_surface=True, nt=nt, B=1,
                     memory=BoundarySaving(storage="gpu"))

torch.cuda.reset_peak_memory_stats()
vp_ref = torch.tensor(vp, device=dev, requires_grad=True)
rec_ref = build(vp.shape)(wav, src, rec, models=[vp_ref])
(0.5 * rec_ref.pow(2).sum()).backward()
g_ref = vp_ref.grad.detach().cpu().numpy()
peak1 = torch.cuda.max_memory_allocated() / 2**30
print(f"single domain: record {tuple(rec_ref.shape)}  peak {peak1:.2f} GiB")
np.save("/tmp/ot_ref_grad.npy", g_ref)

from sweep.equations import Acoustic3D from sweep.propagator.torch import PropTorch from sweep.propagator.options import BoundarySaving # 5 Hz Ricker + order-8: fmax (~12.5 Hz) stays above 3.5 ppw at vmin -> no visible grid dispersion dev, so, abcn, nt, DT = torch.device("cuda:0"), 8, 20, 1500, 0.002 def ricker(nt, dt, fm=5.0, delay=0.25): t = np.arange(nt, dtype=np.float32) * dt - delay a = np.pi * fm * t return ((1.0 - 2.0 * a**2) * np.exp(-(a**2))).astype(np.float32) wav = torch.as_tensor(ricker(nt, DT), device=dev) src = np.array([[[nx//2, ny//2, 4]]], dtype=np.int64) # (x, y, z) rec = np.array([[[ix, ny//2, 4] for ix in range(6, nx-6, 2)]], dtype=np.int64) def build(shape): return PropTorch(Acoustic3D(spatial_order=so, device=dev, backend="torch"), backend="torch", impl="c", shape=shape, dev=dev, dh=dh, dt=DT, source_type=["h1"], receiver_type=["h1"], abcn=abcn, free_surface=True, nt=nt, B=1, memory=BoundarySaving(storage="gpu")) torch.cuda.reset_peak_memory_stats() vp_ref = torch.tensor(vp, device=dev, requires_grad=True) rec_ref = build(vp.shape)(wav, src, rec, models=[vp_ref]) (0.5 * rec_ref.pow(2).sum()).backward() g_ref = vp_ref.grad.detach().cpu().numpy() peak1 = torch.cuda.max_memory_allocated() / 2**30 print(f"single domain: record {tuple(rec_ref.shape)} peak {peak1:.2f} GiB") np.save("/tmp/ot_ref_grad.npy", g_ref)

single domain: record (1, 1500, 196, 1)  peak 16.14 GiB

Multi-tile DD¶

Launched as a subprocess — one rank per GPU. The worker re-runs the same shot under ModelParallel, all-reduces the gradient, and compares against the single-domain arrays saved above.

In [4]:

worker = r"""
import os, sys, numpy as np, torch, torch.distributed as dist
from sweep.datasets import load
from sweep.equations import Acoustic3D
from sweep.parallel import MeshTopology, pad_to_mesh
from sweep.parallel.dd_propagator import ModelParallel
from sweep.propagator.torch import PropTorch
from sweep.propagator.options import BoundarySaving

dist.init_process_group("nccl")
rank, world = dist.get_rank(), dist.get_world_size()
li = int(os.environ.get("LOCAL_RANK", rank)) % max(1, torch.cuda.device_count())
torch.cuda.set_device(li); dev = torch.device(f"cuda:{li}")

py, px = int(sys.argv[1]), int(sys.argv[2])
d  = load("overthrust", variant="3d-acoustic", downsample=2)
vp = np.ascontiguousarray(d["vp"]).astype("float32"); dh = float(d["dh"][0])
TILE = 4                                   # must match the notebook's pad
vp = pad_to_mesh(vp, py=TILE, px=TILE)
nz, ny, nx = vp.shape
so, abcn, nt, DT = 8, 20, 1500, 0.002

t = np.arange(nt, dtype=np.float32) * DT - 0.25
a = np.pi * 5.0 * t
wav = torch.as_tensor(((1 - 2*a**2) * np.exp(-a**2)).astype("float32"), device=dev)
src = np.array([[[nx//2, ny//2, 4]]], dtype=np.int64)
rec = np.array([[[ix, ny//2, 4] for ix in range(6, nx-6, 2)]], dtype=np.int64)

prop = PropTorch(Acoustic3D(spatial_order=so, device=dev, backend="torch"),
                 backend="torch", impl="c", shape=vp.shape, dev=dev, dh=dh, dt=DT,
                 source_type=["h1"], receiver_type=["h1"], abcn=abcn,
                 free_surface=True, nt=nt, B=1,
                 memory=BoundarySaving(storage="gpu"))
mesh = MeshTopology(py=py, px=px, shot_groups=1, world_size=world, rank=rank)
ddp  = ModelParallel(prop, mesh)

torch.cuda.reset_peak_memory_stats()
vp_dd = torch.tensor(vp, device=dev, requires_grad=True)
rec_dd = ddp(wav, src, rec, models=[vp_dd])
(0.5 * rec_dd.pow(2).sum()).backward()
dist.all_reduce(vp_dd.grad)
peak = torch.cuda.max_memory_allocated() / 2**30

if rank == 0:
    g_ref = np.load("/tmp/ot_ref_grad.npy")
    g_dd  = vp_dd.grad.detach().cpu().numpy()
    dif   = np.abs(g_dd - g_ref)
    print(f"DD_RESULT py={py} px={px} peak_per_tile={peak:.2f}GiB "
          f"grad_bitexact={np.array_equal(g_dd, g_ref)} "
          f"max_abs_diff={dif.max():.3e} "
          f"rel_l2={np.linalg.norm(g_dd-g_ref)/(np.linalg.norm(g_ref)+1e-30):.3e}")
    np.save("/tmp/ot_dd_grad.npy", g_dd)
dist.barrier(); dist.destroy_process_group()
"""
import os, sys, subprocess, pathlib
pathlib.Path("/tmp/ot_dd_worker.py").write_text(worker)

ngpu = min(torch.cuda.device_count(), 4)
py, px = (2, 2) if ngpu == 4 else (1, ngpu)          # 4 GPUs -> both a y- and an x-cut
print(f"{ngpu}-way DD  (py={py}, px={px})")
env = dict(os.environ, PYTHONPATH=str(pathlib.Path(sweep.__file__).resolve().parents[1]))
res = subprocess.run([sys.executable, "-m", "torch.distributed.run", "--standalone",
                      f"--nproc-per-node={ngpu}", "/tmp/ot_dd_worker.py", str(py), str(px)],
                     capture_output=True, text=True, env=env)
out = res.stdout + res.stderr
for ln in out.splitlines():
    if "DD_RESULT" in ln:
        print(ln)
if res.returncode != 0:                     # show what actually went wrong
    print("--- subprocess failed, tail of its output ---")
    print("\n".join(out.strip().splitlines()[-15:]))
print("exit:", res.returncode)

worker = r""" import os, sys, numpy as np, torch, torch.distributed as dist from sweep.datasets import load from sweep.equations import Acoustic3D from sweep.parallel import MeshTopology, pad_to_mesh from sweep.parallel.dd_propagator import ModelParallel from sweep.propagator.torch import PropTorch from sweep.propagator.options import BoundarySaving dist.init_process_group("nccl") rank, world = dist.get_rank(), dist.get_world_size() li = int(os.environ.get("LOCAL_RANK", rank)) % max(1, torch.cuda.device_count()) torch.cuda.set_device(li); dev = torch.device(f"cuda:{li}") py, px = int(sys.argv[1]), int(sys.argv[2]) d = load("overthrust", variant="3d-acoustic", downsample=2) vp = np.ascontiguousarray(d["vp"]).astype("float32"); dh = float(d["dh"][0]) TILE = 4 # must match the notebook's pad vp = pad_to_mesh(vp, py=TILE, px=TILE) nz, ny, nx = vp.shape so, abcn, nt, DT = 8, 20, 1500, 0.002 t = np.arange(nt, dtype=np.float32) * DT - 0.25 a = np.pi * 5.0 * t wav = torch.as_tensor(((1 - 2*a**2) * np.exp(-a**2)).astype("float32"), device=dev) src = np.array([[[nx//2, ny//2, 4]]], dtype=np.int64) rec = np.array([[[ix, ny//2, 4] for ix in range(6, nx-6, 2)]], dtype=np.int64) prop = PropTorch(Acoustic3D(spatial_order=so, device=dev, backend="torch"), backend="torch", impl="c", shape=vp.shape, dev=dev, dh=dh, dt=DT, source_type=["h1"], receiver_type=["h1"], abcn=abcn, free_surface=True, nt=nt, B=1, memory=BoundarySaving(storage="gpu")) mesh = MeshTopology(py=py, px=px, shot_groups=1, world_size=world, rank=rank) ddp = ModelParallel(prop, mesh) torch.cuda.reset_peak_memory_stats() vp_dd = torch.tensor(vp, device=dev, requires_grad=True) rec_dd = ddp(wav, src, rec, models=[vp_dd]) (0.5 * rec_dd.pow(2).sum()).backward() dist.all_reduce(vp_dd.grad) peak = torch.cuda.max_memory_allocated() / 2**30 if rank == 0: g_ref = np.load("/tmp/ot_ref_grad.npy") g_dd = vp_dd.grad.detach().cpu().numpy() dif = np.abs(g_dd - g_ref) print(f"DD_RESULT py={py} px={px} peak_per_tile={peak:.2f}GiB " f"grad_bitexact={np.array_equal(g_dd, g_ref)} " f"max_abs_diff={dif.max():.3e} " f"rel_l2={np.linalg.norm(g_dd-g_ref)/(np.linalg.norm(g_ref)+1e-30):.3e}") np.save("/tmp/ot_dd_grad.npy", g_dd) dist.barrier(); dist.destroy_process_group() """ import os, sys, subprocess, pathlib pathlib.Path("/tmp/ot_dd_worker.py").write_text(worker) ngpu = min(torch.cuda.device_count(), 4) py, px = (2, 2) if ngpu == 4 else (1, ngpu) # 4 GPUs -> both a y- and an x-cut print(f"{ngpu}-way DD (py={py}, px={px})") env = dict(os.environ, PYTHONPATH=str(pathlib.Path(sweep.__file__).resolve().parents[1])) res = subprocess.run([sys.executable, "-m", "torch.distributed.run", "--standalone", f"--nproc-per-node={ngpu}", "/tmp/ot_dd_worker.py", str(py), str(px)], capture_output=True, text=True, env=env) out = res.stdout + res.stderr for ln in out.splitlines(): if "DD_RESULT" in ln: print(ln) if res.returncode != 0: # show what actually went wrong print("--- subprocess failed, tail of its output ---") print("\n".join(out.strip().splitlines()[-15:])) print("exit:", res.returncode)

4-way DD  (py=2, px=2)

DD_RESULT py=2 px=2 peak_per_tile=4.29GiB grad_bitexact=True max_abs_diff=0.000e+00 rel_l2=0.000e+00
exit: 0

Does the gradient survive the cuts?¶

A halo or PML bug at a tile boundary shows up as a stripe on the cut, which a global norm would average away. Slice all three ways — depth, inline and crossline, straight across the cut planes — and plot the difference.

In [5]:

from matplotlib.colors import TwoSlopeNorm

g_dd = np.load("/tmp/ot_dd_grad.npy")
cx, cy = nx // px, ny // py
cuts_x = [i * cx for i in range(1, px)]
cuts_y = [i * cy for i in range(1, py)]
iz, iy, ix = nz // 3, ny // 2, nx // 2      # iy/ix pass through the source

rows = [  # (slice, title, xlabel, ylabel, vlines, hlines, clip_from)
    (lambda g: g[iz],       f"depth z={iz}",                 "x", "y", cuts_x, cuts_y, 0),
    (lambda g: g[:, iy, :], f"inline y={iy} (clip z>12)",    "x", "z", cuts_x, [],     12),
    (lambda g: g[:, :, ix], f"crossline x={ix} (clip z>12)", "y", "z", cuts_y, [],     12),
]
fig, ax = plt.subplots(3, 3, figsize=(15, 10), constrained_layout=True,
                       gridspec_kw={"height_ratios": [2.2, 1, 1]})
for r_, (sl, ttl, xl, yl, vlines, hlines, cf) in enumerate(rows):
    a, b = sl(g_ref), sl(g_dd)
    d = a[cf:]
    lo, hi = np.percentile(d, [0.5, 99.5])   # asymmetric: follow the real range
    m = np.abs(d).max() or 1.0               # keep 0 strictly inside the range
    nrm = TwoSlopeNorm(vmin=min(lo, -1e-3 * m), vcenter=0.0, vmax=max(hi, 1e-3 * m))
    for k, (img, name) in enumerate(((a, "single"), (b, "DD"), (b - a, "DD - single"))):
        im = ax[r_, k].imshow(img, cmap="seismic", norm=nrm,
                              aspect="auto", origin="upper")
        for c in vlines: ax[r_, k].axvline(c, color="k", ls="--", lw=1.0)
        for c in hlines: ax[r_, k].axhline(c, color="k", ls="--", lw=1.0)
        ax[r_, k].set_title(f"{name}   {ttl}", fontsize=10)
        ax[r_, k].set_xlabel(xl); ax[r_, k].set_ylabel(yl)
        plt.colorbar(im, ax=ax[r_, k], shrink=0.85)
plt.show()

iy_t = cuts_y[0] if cuts_y else ny // 2
fig, ax = plt.subplots(2, 1, figsize=(11, 5.6), sharex=True, constrained_layout=True)
ax[0].plot(g_ref[iz, iy_t], lw=2.4, color="0.55", label="single")
ax[0].plot(g_dd[iz, iy_t], lw=0.9, color="crimson", ls="--", label="DD")
ax[0].set_ylabel("grad"); ax[0].legend(fontsize=9)
ax[0].set_title(f"trace across the cuts   z={iz}, y={iy_t}", fontsize=10)
r = g_dd[iz, iy_t] - g_ref[iz, iy_t]
ax[1].plot(r, lw=1.0, color="tab:blue")
ax[1].set_xlabel("x"); ax[1].set_ylabel("DD - single")
ax[1].set_title(f"residual, max={np.abs(r).max():.3e}", fontsize=10)
if np.abs(r).max() == 0.0: ax[1].set_ylim(-1, 1)
for a_ in ax:
    for c in cuts_x: a_.axvline(c, color="k", ls=":", lw=1.2)
plt.show()

print(f"grad bit-exact : {np.array_equal(g_dd, g_ref)}")
print(f"max |DD-single|: {np.abs(g_dd-g_ref).max():.3e}")

from matplotlib.colors import TwoSlopeNorm g_dd = np.load("/tmp/ot_dd_grad.npy") cx, cy = nx // px, ny // py cuts_x = [i * cx for i in range(1, px)] cuts_y = [i * cy for i in range(1, py)] iz, iy, ix = nz // 3, ny // 2, nx // 2 # iy/ix pass through the source rows = [ # (slice, title, xlabel, ylabel, vlines, hlines, clip_from) (lambda g: g[iz], f"depth z={iz}", "x", "y", cuts_x, cuts_y, 0), (lambda g: g[:, iy, :], f"inline y={iy} (clip z>12)", "x", "z", cuts_x, [], 12), (lambda g: g[:, :, ix], f"crossline x={ix} (clip z>12)", "y", "z", cuts_y, [], 12), ] fig, ax = plt.subplots(3, 3, figsize=(15, 10), constrained_layout=True, gridspec_kw={"height_ratios": [2.2, 1, 1]}) for r_, (sl, ttl, xl, yl, vlines, hlines, cf) in enumerate(rows): a, b = sl(g_ref), sl(g_dd) d = a[cf:] lo, hi = np.percentile(d, [0.5, 99.5]) # asymmetric: follow the real range m = np.abs(d).max() or 1.0 # keep 0 strictly inside the range nrm = TwoSlopeNorm(vmin=min(lo, -1e-3 * m), vcenter=0.0, vmax=max(hi, 1e-3 * m)) for k, (img, name) in enumerate(((a, "single"), (b, "DD"), (b - a, "DD - single"))): im = ax[r_, k].imshow(img, cmap="seismic", norm=nrm, aspect="auto", origin="upper") for c in vlines: ax[r_, k].axvline(c, color="k", ls="--", lw=1.0) for c in hlines: ax[r_, k].axhline(c, color="k", ls="--", lw=1.0) ax[r_, k].set_title(f"{name} {ttl}", fontsize=10) ax[r_, k].set_xlabel(xl); ax[r_, k].set_ylabel(yl) plt.colorbar(im, ax=ax[r_, k], shrink=0.85) plt.show() iy_t = cuts_y[0] if cuts_y else ny // 2 fig, ax = plt.subplots(2, 1, figsize=(11, 5.6), sharex=True, constrained_layout=True) ax[0].plot(g_ref[iz, iy_t], lw=2.4, color="0.55", label="single") ax[0].plot(g_dd[iz, iy_t], lw=0.9, color="crimson", ls="--", label="DD") ax[0].set_ylabel("grad"); ax[0].legend(fontsize=9) ax[0].set_title(f"trace across the cuts z={iz}, y={iy_t}", fontsize=10) r = g_dd[iz, iy_t] - g_ref[iz, iy_t] ax[1].plot(r, lw=1.0, color="tab:blue") ax[1].set_xlabel("x"); ax[1].set_ylabel("DD - single") ax[1].set_title(f"residual, max={np.abs(r).max():.3e}", fontsize=10) if np.abs(r).max() == 0.0: ax[1].set_ylim(-1, 1) for a_ in ax: for c in cuts_x: a_.axvline(c, color="k", ls=":", lw=1.2) plt.show() print(f"grad bit-exact : {np.array_equal(g_dd, g_ref)}") print(f"max |DD-single|: {np.abs(g_dd-g_ref).max():.3e}")

grad bit-exact : True
max |DD-single|: 0.000e+00

What DD bought¶

peak_per_tile above is the per-GPU peak; compare it with the single-domain peak printed earlier. DD trades a halo exchange per step for a model (and wavefield, and boundary ring) that no longer has to fit on one card — which is what makes 3-D field-scale grids reachable.

Two limits worth remembering:

• Comm/compute overlap is x-cut only. With a y-cut in the mesh the phased forward is switched off automatically and the serial path is used.
• A source within M cells of an x-cut also disables overlap (it would be injected after the strips were exchanged). Correctness is preserved either way — see test/dd_src_on_cut_check.py.