sp-DBA: a general framework for adaptive transform-domain computation

sp-DBA: a general framework for adaptive transform-domain computation Jingkun Jiang

Pingchuan Deng

Yang Xia*

College of Materials Science and Engineering, Hunan University Changsha 410082, Hunan, People’s Republic of China

arXiv:2609.03922v1 [cs.CE] 3 Sep 2026

Email: JJJK0506@hnu.edu.cn

*

Corresponding author: yxia@hnu.edu.cn

Abstract Transform-domain methods simplify analysis and computation, making them central to scientific computing and signal processing. However, existing adaptive strategies often introduce new data structures or require substantial workflow redesign, limiting efficient execution on massively parallel hardware. Here we present spectral dynamic block activation (sp-DBA), an adaptive acceleration framework for transform-domain workflows. sp-DBA dynamically activates transform-domain computation blocks as a calculation evolves, concentrating computation where it is needed while retaining numerical accuracy and optimized transform operations. Across representative workflows in materials science, biology, and optical communications, the demonstrated implementations achieve speedups of up to 28.1-fold for transform-domain updates and up to 8.4-fold overall acceleration; strongand weak-scaling tests on up to eight GPUs show that the local adaptive update remains effective within distributed FFT workflows. By introducing execution-level adaptivity into mature transform-domain workflows without rebuilding existing solvers, sp-DBA extends the scale, duration, and complexity of scientific simulations and signal-processing calculations on modern parallel systems.

1

Contents 1 Introduction

3

2 Methods

4

3 Results 6 3.1 Polycrystalline coarsening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.2 Multicomponent pattern formation . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.3 Digital backpropagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.4 Performance and scalability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4 Discussion

13

References

14

Supplementary Information

20

S1 sp-DBA framework S1.1 Boundary-aware transition smoothing . . . . . . . . . . . . . . . . . . . . . . . . S1.2 Neighbor activation in half-spectrum FFT layouts . . . . . . . . . . . . . . . . . S1.3 Dynamic threshold control and periodic refresh . . . . . . . . . . . . . . . . . . .

21 21 23 23

S2 sp-DBA integration with transform-domain solvers S2.1 Integration procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S2.2 Implementation overhead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S2.3 Comparison with adaptive spectral and wavelet methods . . . . . . . . . . . . . .

25 26 27 27

S3 XPFC polycrystalline coarsening S3.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S3.2 Numerical scheme and simulation configuration . . . . . . . . . . . . . . . . . . . S3.3 Grain reconstruction and system-size comparison . . . . . . . . . . . . . . . . . .

28 28 29 31

S4 Reaction–cross-diffusion pattern formation

32

S5 Digital backpropagation

35

S6 Numerical validation

38

S7 Single-GPU performance evaluation 40 S7.1 PFC cost decomposition and acceleration . . . . . . . . . . . . . . . . . . . . . . 41 S7.2 DBP cost decomposition and acceleration . . . . . . . . . . . . . . . . . . . . . . 42 S8 Multi-GPU sp-DBA with distributed FFTs S8.1 Distributed FFT workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S8.2 Strong and weak scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S8.3 CUDA-aware MPI communication . . . . . . . . . . . . . . . . . . . . . . . . . . S8.4 Activation stability under distributed execution . . . . . . . . . . . . . . . . . . .

42 43 45 46 46

S9 Visualization

47

2

1 Introduction Transform-domain methods are central to scientific computing, engineering, and signal processing because they express fields and signals in representations better suited to analysis and numerical computation [1, 2]. Among these methods, Fourier-based workflows are especially widespread. In a Fourier representation, numerical operations such as differentiation and convolution take simpler forms [3], while fast Fourier transform (FFT) algorithms substantially reduce the cost of repeated transforms [4]. With advances in parallel computing hardware, highly optimized FFT libraries have been developed for CPUs [5] and GPUs [6], while distributed FFT implementations have enabled large-scale parallel transforms [7, 8, 9]. As a result, Fourier-based workflows are now widely used in areas including materials simulation [10, 11], fluid simulation [12], computational physics [13, 14], medical imaging [15], optical signal processing [16, 17], and operator learning [18, 19]. In addition to advances in transform algorithms and implementations, the increasing computational demands of large-scale and complex scientific problems call for adaptive methods that respond to the evolving structure of fields and signals while avoiding unnecessary computation. In transform-domain computation, much of the prior work realizes this adaptivity by changing the numerical representation of the field or signal. Examples include sparse Fourier methods [20, 21], wavelet and multiresolution methods [22, 23, 24], and adaptive spectral or reduced-order approaches [25, 26]. However, in mature Fourier-based workflows built around regular arrays and highly optimized FFT libraries, such changes often require reorganizing the solution data or introducing specialized data structures for adaptivity, limiting direct reuse of existing implementations [5, 7, 27]. Moreover, this difficulty is amplified on GPUs and other parallel hardware, where irregular execution and data access often reduce the efficiency of adaptive workloads [28, 29, 30]. For efficient execution on massively parallel architectures, adaptive methods commonly organize computation into patches or blocks that serve as regular units of parallel work [31, 32, 33]. A recent practical realization of this strategy is dynamic block activation (DBA), which organizes computation into blocks and dynamically activates them according to evolving features of the physical model [34]. By preserving a regular blockwise execution pattern on GPUs, DBA offers a straightforward approach to accelerating existing continuum models through adaptivity. However, applying blockwise activation directly to transform-domain workflows presents two fundamental difficulties. First, before transformation, the data are represented by local field values, whereas after transformation, they are represented by modes whose contributions span the physical domain, so a transformed-field block does not correspond directly to a local region in the original field [1, 2, 3]. Second, a typical transform-domain workflow tightly couples forward and inverse transforms, the transformed-field update, and other numerical operations, so adaptive decisions must remain consistent with the full workflow to preserve numerical stability and accuracy [2, 3]. Consequently, the challenge is to introduce adaptivity without sacrificing the numerical structure and efficient parallel implementations that make mature transform-domain workflows practical. In this work, we present spectral dynamic block activation (sp-DBA), an adaptive acceleration framework for transform-domain workflows. sp-DBA applies blockwise activation to transformedfield computation, coordinating adaptive execution with the surrounding transform workflow on parallel hardware. With only minor changes, we demonstrate the same framework in three representative applications, including polycrystalline coarsening, multicomponent reaction– cross-diffusion, and digital backpropagation, spanning materials science, biology, and optical communications. Across the three implementations, numerical accuracy is maintained while transform-domain update speedups exceed twentyfold in the strongest configuration. We further evaluate stage-level cost decomposition and distributed PFC execution on 1–8 GPUs, including strong and weak scaling, to identify where local acceleration persists and where global FFT communication limits end-to-end gains (Fig. 5; Supplementary Sections S7 and S8). More importantly, this acceleration expands the scientific questions accessible within a fixed computing 3

budget. Larger domains and longer evolutions provide stronger statistics and clearer separation of intrinsic behavior from finite-size and transient effects, while the same computational savings make richer coupled models and more demanding signal-processing studies accessible.

2 Methods

a

Fourier space

Amplitude

ky

kx

Data layer

Block state

Transition

Active

Frozen

Algorithm layer

GPU execution

Hardware execution

1 sp-DBA block = 1 thread block

Shared memory states neighbors

Device memory

b

c

Transform workflow and sp-DBA Pre-transform

State preparation

GPU execution pattern Pointwise

Blockwise

Threads

Threads

Device memory

Device memory

Threshold control

Transform

State reading

Blockwise update Transformed-field

update update core Inverse transform

Block activation Apply update core A

T

F

Consistency control Post-transform

Consistency

Refresh

Figure 1: Overview of sp-DBA in transform-domain solver workflows. a, Data, algorithm, and hardware mapping. The transformed field is partitioned into active, transition, and frozen blocks for full updates, smoothed updates, and update skipping, respectively. Each sp-DBA block maps to one CUDA thread block, with block-state information stored in shared memory. b, Transform workflow and sp-DBA insertion. sp-DBA contains three modules: state preparation, blockwise update, and consistency control, which respectively manage block states, apply updates adaptively, and maintain numerical consistency. c, GPU execution patterns. Pointwise activation assigns different update paths to threads within one CUDA thread block, whereas blockwise activation assigns one path to the full thread block.

4

Spectral dynamic block activation. sp-DBA is an adaptive acceleration framework for transform-domain workflows. In the implementations presented here, GPUs serve as a representative parallel platform; their block-based execution model provides a natural mapping for sp-DBA. Throughout this work, Regular denotes the full-update mode with sp-DBA disabled, in which the model-specific update core is applied to the full transformed field. With sp-DBA enabled, the same update core is applied according to block state, while the existing transform workflow remains in place. The framework consists of three modules: state preparation, blockwise update, and consistency control. For blockwise execution, the transformed-field data are partitioned into sp-DBA blocks, each serving as an activation unit and mapping one-to-one to a CUDA thread block (Fig. 1a). State preparation evaluates an application-specific activation measure for each block and compares it with a threshold. The threshold may be fixed, as in DBA [34], or adjusted dynamically to maintain a target average active-block ratio, which includes both active and transition blocks. Blocks meeting the threshold condition are marked active. Among the remaining blocks, those neighboring an active block enter the transition state, whereas the rest are frozen. The block state then determines the update: active blocks execute the model-specific update core, transition blocks apply a boundary-aware weighted update, and frozen blocks skip the update and retain their stored transformed-field values until reactivation or periodic refresh (Supplementary Section S1.3). During GPU execution, the current and neighboring block states, together with the smoothing table, are loaded into shared memory to evaluate the transition weight locally (Supplementary Section S1.1). A representative step in a transform-domain workflow includes a pre-transform stage, a forward transform, a transformed-field update, an inverse transform, and a post-transform stage (Fig. 1b). For example, in a pseudo-spectral solver, the pre-transform stage evaluates nonlinear or reaction terms in real space, and the forward transform maps them to Fourier space. At the transformed-field update stage, sp-DBA performs state preparation, blockwise update, and consistency control. The inverse transform and post-transform stages then proceed as in Regular. At the hardware level, blockwise activation assigns one update decision to each CUDA thread block, whereas pointwise activation assigns decisions to individual transformed-field values (Fig. 1c). These design choices provide four practical benefits: hardware-aligned execution, complementary acceleration, solver compatibility, and simple integration. 1. Acceleration with hardware-aligned block execution. Assigning one update decision to each CUDA thread block lets frozen blocks avoid unnecessary updates while keeping all threads in the block on a common execution path (Fig. 1c). 2. Complementary acceleration with optimized transforms. Optimized transform implementations, including FFTW [5], cuFFT [6], heFFTe [8], and cuFFTMp [9], reduce the cost of the forward and inverse transforms. sp-DBA complements these optimized transforms and reduces the number of blocks executing the model-specific update, allowing both forms of acceleration to be combined within the same workflow. 3. Compatibility with different numerical solvers. sp-DBA operates on each solver’s existing transformed-field update while leaving its governing equations, numerical discretization, transformed representation, and transform workflow unchanged. The same framework is demonstrated in PFC, RCD, and DBP despite differences in their equations and transformedfield layouts (Supplementary Section S2.1). By contrast, representative adaptive spectral and wavelet methods change the basis or computational grid and therefore require state transfer or dedicated hierarchical data structures [25, 27] (Supplementary Section S2.3). 4. Simple integration into existing solvers. Common block operations are implemented once in the shared framework, while each solver reuses its existing update core and supplies the 5

required model-specific data and control settings [35]. Direct integration into the three solvers adds code amounting to about 8–16% of the corresponding Regular solver core. Including model-adapter code raises the total to about 13–28%; the shared framework is excluded from both per-solver counts (Supplementary Section S2.2).

3 Results 3.1 Polycrystalline coarsening a

t Regular

sp-DBA

sp-DBA

Atomic scale structure

Real space Fourier space

t = 10

b

Amplitude

sp-DBA state

0°

4.5e-2

Active blocks Transition blocks

4.5e-4 4.5e-6

t = 100

c

60°

Orientation θ t = 10,000

Grain boundary/defects

d

Figure 2: Polycrystalline coarsening with sp-DBA. a, Layered evolution of XPFC polycrystalline coarsening. From left to right, the columns show snapshots at dimensionless times t = 10, 100, and 10, 000, with Regular used for the first and sp-DBA for the latter two, all on a 20482 grid. From top to bottom, the rows show real-space morphology colored by local crystal orientation θ, a three-dimensional surface plot of the Fourier-space amplitude with the zero-frequency component omitted, and sp-DBA block states. Red and yellow mark active and transition blocks, respectively, and the inset highlights grain boundaries and atomic-scale defects. b, Normalized grain-size distributions for Regular and sp-DBA during late-stage coarsening on an 81922 grid. c, Free energy F [ψ] and mean density ψ̄ for Regular and sp-DBA during coarsening on an 81922 grid. d, Fourier-space update speedup relative to Regular as a function of grid size. Percentages indicate target average active-block ratios maintained by dynamic threshold control, and dashed lines show the corresponding ideal speedups.

We first test sp-DBA in XPFC polycrystalline coarsening, a representative transform-domain workflow linking a lattice-periodic density field to long-time grain-boundary migration. Phasefield-crystal (PFC) models describe crystalline microstructure through a continuous density field ψ(r) with periodic order at atomic length scales, capturing lattice ordering, elasticity, defects, and grain boundaries while accessing diffusive time scales beyond direct atomistic simulations [36, 37]. XPFC extends this framework with a two-point correlation kernel C2 that selects the crystalline length scale and lattice symmetry [38, 39]. As PFC methods have developed, XPFC and related formulations have been applied to crystalline coarsening, structural transformations, and defect

6

evolution [40, 41]. For the XPFC calculations considered here, we use the free energy Z

F [ψ] =

ω u 1 1 dr ψ 2 − ψ 3 + ψ 4 − 2 6 12 2 



ZZ

dr dr ′ ψ(r)C2 (|r − r ′ |)ψ(r ′ ).

(1)

In the pseudo-spectral XPFC solver, the local polynomial terms are evaluated from ψ(r) in real space, whereas the correlation term is evaluated in Fourier space as Ĉ2 (k)ψ̂(k), where Ĉ2 (k) is the Fourier representation of C2 [42]. During coarsening, density peaks, defects, and grain boundaries are distributed throughout the physical domain, whereas the dominant Fourier amplitudes concentrate near the Bragg rings associated with the crystalline length scale (Fig. 2a). Accordingly, sp-DBA reduces or skips updates in blocks away from these rings. This contrast between distributed real-space morphology and localized Fourier amplitudes makes XPFC particularly suitable for testing sp-DBA. Polycrystalline states are initialized from randomly oriented crystalline seeds. All Fourier blocks are kept active during the rapidly evolving early growth stage; once the polycrystalline structure forms and grain-boundary migration becomes dominant, sp-DBA is enabled for the subsequent stage shown in Fig. 2a. For each Fourier block Bj , the block energy EBj = P 2 k∈Bj |ψ̂(k)| identifies the blocks requiring updates near the Bragg rings. Their gradual evolution allows dynamic threshold control to maintain target average active-block ratios of 5%, 10%, 15%, and 20% for controlled performance comparisons (Supplementary Section S3). To verify agreement with Regular, we compare grain statistics and thermodynamic quantities. Across the tested target ratios, the normalized grain-size distribution (GSD) [43] closely follows the Regular result (Fig. 2b), consistent with grain-growth analyses based on equivalent grain diameter [44, 45]. In addition, the free energy F [ψ] decreases similarly in the Regular and sp-DBA calculations, while the mean density ψ̄ is conserved (Fig. 2c). Agreement across grain statistics, free-energy evolution, and mean density confirms that sp-DBA preserves the coarsening statistics and thermodynamic behavior under the tested activation settings. For performance, lower target active-block ratios produce larger Fourier-space update speedups. For N ≥ 2048, the 20%, 10%, and 5% targets yield approximately 5×, 10×, and 18–20× update speedups, respectively (Fig. 2d). On smaller grids, by contrast, the Bragg-ring region occupies a larger fraction of the Fourier blocks, leaving fewer blocks to skip, while block-state operations account for a larger share of the update time. With increasing grid size, the measured speedups approach the corresponding ideal values, indicating low block-state overhead. This growing benefit at large domains is particularly relevant to coarsening studies, where retaining sufficiently large grain populations is essential for robust late-stage statistics.

7

3.2 Multicomponent pattern formation

a step = 1,000

step = 4,000

sp-DBA

b

z

Fourier space

step = 20,000

Regular

c

d

~3% active blocks

Active blocks y

x

Transition blocks

Figure 3: Reaction–cross-diffusion pattern formation with sp-DBA. a, Real-space u1 morphology on a 2563 grid with M = 8 at 1,000, 4,000, and 20,000 steps. Blue surfaces show sp-DBA; gray wireframes at 4,000 and 20,000 steps show the Regular mode. b, Fourier-space block states for the same sp-DBA calculation. The average active-block ratio is approximately 3%, counting active and transition blocks. Red and light blue mark active and transition blocks, respectively. c, Radially averaged structure factor S(k) of u1 for Regular and sp-DBA. d, Normalized Fourier-space update cost as the number of components M increases, relative to Regular. Labels give the update-stage speedups of sp-DBA with eigendecomposition.

We next test sp-DBA in reaction–cross-diffusion (RCD) pattern formation, a multicomponent transform-domain workflow in which each stored Fourier mode carries a coupled update involving M components. In classical reaction–diffusion systems, reaction kinetics and diffusion can destabilize a homogeneous state and generate spatial patterns from small perturbations [46]. Cross diffusion extends this setting by allowing the flux of one component to depend on gradients of other components, providing a mechanism for pattern formation in multicomponent systems [47, 48]. More generally, recent analyses show that non-diagonal diffusion matrices can support diffusion-driven instabilities in reaction–cross-diffusion systems [49]. For the calculations considered here, we use an RCD system in which a constant cross-diffusion matrix couples the M components, ∂u = ∇ · (D∇u) + R(u), (2) ∂t where u = (u1 , . . . , uM )T , D is a constant cross-diffusion matrix, and R(u) contains the nonlinear reaction terms (Supplementary Section S4). In the Fourier implementation used here, R(u) is evaluated in real space, whereas the crossdiffusion term is updated in Fourier space [1, 2]. At each Fourier mode k, the transformed 8

state is the component vector û(k), and the matrix D couples its M entries. Thus, each stored Fourier mode carries a coupled M × M update rather than a scalar multiplication. The eigendecomposition of the constant matrix D is computed once before time stepping and reused during the Fourier update [50] (Supplementary Section S4.3). For the RCD performance comparisons, Regular and sp-DBA use the direct coupled update unless eigendecomposition is stated explicitly. RCD simulations are initialized with small random perturbations around a homogeneous reference state, after which unstable modes grow and saturate into a labyrinthine pattern. For visualization and spectral comparison, we use component u1 for the morphology and structure factor S(k). For block activation, sp-DBA combines all components through the root-meanP 2 1/2 . Accordingly, the same block state is square Fourier amplitude A(k) = (M −1 M i=1 |ûi (k)| ) applied to all M components; active and transition blocks update the full component vector, whereas frozen blocks skip the coupled operation (Supplementary Section S4.2). To assess agreement with Regular, we compare real-space morphology and wavelength selection. As the instability grows and saturates, the u1 (r) morphology from sp-DBA closely matches the Regular result, with the blue sp-DBA surfaces overlapping the gray Regular wireframes at 4,000 and 20,000 steps (Fig. 3a). In Fourier space, active and transition blocks concentrate in the low-wavenumber band where the unstable modes grow, leaving most stored Fourier blocks frozen (Fig. 3b). This localization is consistent with the finite instability band that selects the dominant wavelength [49, 51]. To quantify wavelength selection, we compute the radially averaged structure factor S(k) = ⟨S(k)⟩|k|=k from u1 . The peak positions and amplitudes agree closely between Regular and sp-DBA (Fig. 3c). Outside the unstable band, high-wavenumber modes decay rapidly, so the required updates are concentrated at low wavenumbers (Supplementary Section S4.2). The agreement in morphology and S(k) shows that sp-DBA preserves wavelength selection even though most high-wavenumber blocks are frozen. For performance, increasing M makes each frozen block skip a larger coupled Fourier-space update. Regular applies the M × M operation at every stored wavenumber, whereas spDBA executes it only in active and transition blocks. Relative to Regular, sp-DBA with eigendecomposition gives update-stage speedups of 9.6× at M = 24, 28.1× at M = 48, and 26.7× at M = 64 (Fig. 3d). At smaller M , block-state reading and updates account for a larger fraction of the update-stage cost, resulting in smaller measured speedups. These measurements show that block activation removes a large fraction of the coupled update cost, directly targeting the growing computational burden of high-component models. Fig. 5d and Supplementary Section S4.3 separate this contribution from eigendecomposition.

9

3.3 Digital backpropagation

a

Input transmission Light signal

Fiber

Recovered

DBP split-step loop

Frequency

Frequency

Time

Reference

Fiber output (impaired signal)

sp-DBA

Mode coupling Dispersion Nonlinearity

b

Recovered output

Signal recovery with sp-DBA

sp-DBA

Retained

Spectrum evolution

FFT/IFFT

e litud

mp al a ign

S Frequency index

Recovery distance Amplitude (a.u.) 0.1

c

1

10

100

bloc

-k 1 k st

1000

d

0 ate

k1

v

co

Re

ery

ce

tan

dis

Active Transition

e

Figure 4: Digital backpropagation with sp-DBA. a, Split-step digital backpropagation for optical signal recovery, with sp-DBA applied to the frequency-domain linear update. b, Evolution of the signal spectrum and block states, together with a schematic view of waveform recovery. Red and yellow mark active and transition blocks, respectively. c, Field correlation coefficient of the recovered signal relative to the transmitted reference versus signal length Nt . d, Normalized mean-squared error (NMSE) of the recovered signal versus Nt . In c,d, percentages indicate target average active-block ratios, and the dashed line in d marks the −20 dB reference. e, Normalized frequency-domain update cost for Regular and sp-DBA versus Nt . The default sp-DBA setting targets a 10% average active-block ratio, and labels give the update-stage speedups.

To test sp-DBA in an engineering workflow dominated by repeated frequency-domain updates, we apply it to space-division-multiplexed digital backpropagation (SDM-DBP) for optical signal recovery (Fig. 4a). In SDM links, signals propagate through multiple spatial modes or fiber cores, increasing transmission capacity while producing coupled complex fields affected by dispersion, mode coupling, attenuation, and Kerr nonlinearity [52, 53, 54]. At the receiver, DBP compensates these impairments by numerically reversing signal propagation through split-step updates that alternate frequency-domain linear propagation and time-domain nonlinear correction. Within this split-step procedure, increasing the number of propagation steps generally improves the approximation of inverse propagation, but each added step repeats both updates; longer signal 10

windows and more spatial channels further increase the computational workload [16, 17, 55]. Accordingly, the repeated split-step workload has motivated reduced-complexity methods in optical signal processing [56, 57]. For the coupled field envelope E = [E1 , . . . , EM ]T , the SDM-DBP model used here is ∂E = (D̂ + N̂ )E = ∂z

"

α β2 ∂ 2 − I −i I + iK 2 2 ∂t2

!

8 + iγ 9

M X

! # 2

|Em |

I E,

(3)

m=1

where the linear operator D̂ includes attenuation α, chromatic dispersion β2 , and mode coupling K, while the nonlinear operator N̂ represents the Kerr response with coefficient γ driven by the total modal power (Supplementary Section S5.1). In the numerical implementation, sp-DBA is applied to the repeated frequency-domain evaluation of D̂, while N̂ is evaluated in the time domain [58, 59, 60]. The input field consists of compact Gaussian wave packets with different carrier-frequency offsets, producing spectral bands that remain localized during backpropagation (Fig. 4b). For block activation, sp-DBA P 2 uses the total modal spectral power P(ω) = M m=1 |Ẽm (ω)| . The resulting frequency-block states therefore follow the spectral bands occupied by all modes (Supplementary Section S5.3). To assess recovery quality, we compare each recovered field with the transmitted reference using the field correlation coefficient and normalized mean-squared error (NMSE) over the full multimode complex field. NMSE is commonly used to evaluate recovery in optical DBP and nonlinear-compensation studies [16, 55] (Supplementary Section S5.3). Across the tested signal lengths Nt , the sp-DBA field correlation coefficient remains close to unity, with the lowest plotted value around 0.998 and the larger-Nt cases approaching the Regular result (Fig. 4c). Likewise, the sp-DBA NMSE remains at or below the −20 dB reference, corresponding to a linear NMSE no larger than 10−2 , whereas the Regular calculation gives near-machine-precision recovery for this controlled input (Fig. 4d). Overall, the correlation and NMSE results show that sp-DBA preserves signal-recovery accuracy under the tested activation settings. For performance, increasing Nt expands the sampled frequency grid around these occupied bands, leaving more off-band blocks frozen during the linear update. Consequently, at the default 10% target average active-block ratio, the frequency-domain update-stage speedup increases from 2.2× at Nt = 213 to 21.3× at Nt = 218 (Fig. 4e). Reduced-complexity DBP methods commonly lower cost by narrowing the processed spectrum through sub-bands or bandwidth restrictions [57, 61], simplifying the propagation model [56], or changing the split-step schedule [60, 62]. As the sampled spectrum expands beyond the occupied bands, sp-DBA becomes increasingly effective for longer signal windows, where repeated frequency-domain updates otherwise impose a growing computational cost.

11

3.4 Performance and scalability

a

b

c

d

e

f Regular update

Update time (s)

2

1 9.7%

0

1

sp-DBA update

9.5%

9.0%

8.2%

2

4

8

Number of GPUs

Figure 5: Performance and scalability of sp-DBA. Unless noted otherwise, dynamic threshold control targets a 10% average active-block ratio, and costs in a–c are normalized to the corresponding Regular baseline. FFT, NL, and Upd. denote transform operations, model-specific nonlinear calculations, and transformed-field updates. a, PFC cost decomposition as grid size N increases. b, DBP cost decomposition as signal length Nt increases at M = 32. c, RCD cost decomposition as component number M increases, comparing Regular and sp-DBA with eigendecomposition. d, RCD update time for Regular and sp-DBA with and without eigendecomposition. The inset shows the Regular runtime fractions of the update and all other stages. e, Block-activation estimate and measured combined RCD end-to-end acceleration as a function of M , both relative to Regular. f, PFC update-stage weak scaling at approximately fixed per-GPU workload. Bars show means over three runs, error bars indicate ±1 standard deviation, and labels give measured average active-block ratios.

With transform operations and other model steps retained, end-to-end acceleration depends on the share of the Regular runtime spent in the transformed-field update (Fig. 5a–c). In PFC, the update-stage speedup approaches 10× at large grids, although the two-dimensional FFT and real-space stages still account for a substantial fraction of the runtime (Fig. 5a). In DBP, the resulting update reduction allows end-to-end acceleration to reach 8.4× (Fig. 5b). In RCD, the coupled update accounts for a larger share of the Regular runtime as M increases. Accordingly, the end-to-end acceleration of sp-DBA with eigendecomposition relative to Regular shows an overall upward trend (Fig. 5c). RCD also demonstrates the complementary roles of sp-DBA and eigendecomposition. Eigendecomposition lowers the cost of the matrix operation at each updated Fourier mode, whereas sp-DBA reduces the number of blocks that perform it (Supplementary Section S4.3). At M = 48, combining the two gives the lowest update time and a 28.1× update-stage speedup relative to Regular (Figs. 3d and 5d). For the complete calculation, we compare the measured end-to-end acceleration of the combined configuration with a block-activation estimate (Fig. 5e). The estimate uses the Regular update-stage cost share and average active-block ratio and excludes the eigendecomposition gain. At small M , block-control overhead is more visible, and the measured acceleration is close to or below the estimate. At larger M , eigendecomposition provides a larger additional reduction, and the measured curve rises above the estimate. Finally, we evaluate the sp-DBA update stage within a distributed PFC FFT workflow. 12

Under weak scaling, the per-GPU workload remains approximately fixed, the measured average active-block ratio stays near the 10% target, and the update times of both Regular and sp-DBA remain stable from one to eight GPUs. sp-DBA also reduces the per-GPU update time by about 2–3× (Fig. 5f). Supplementary strong-scaling tests and distributed cost decompositions likewise show lower sp-DBA update times (Supplementary Section S8). The global transpose and all-to-all communication remain unchanged, but the update-stage reduction persists across the tested distributed configurations.

4 Discussion sp-DBA differs from established forms of transform-domain adaptivity in what it adapts. Existing approaches commonly obtain adaptivity by changing the active spectral set [20, 21], the expansion order or basis parameters [25, 63], or the computational grid and coefficient hierarchy [22, 24, 27]. sp-DBA instead introduces adaptivity at the level of execution, leaving the transformed representation and transform workflow unchanged while dynamically controlling where computation is performed. This separation shifts adaptivity from solver-specific numerical redesign to a reusable execution strategy for transform-domain computation. The three demonstrations show that this execution-level principle carries across substantially different scientific and engineering problems. It is notable that the acceleration provided by sp-DBA increases both the spatial scale and evolution time accessible within the same computing budget. For grain-coarsening studies, this broader computational reach is scientifically important because reliable statistics require domains large enough to retain sufficient grains and simulations long enough to distinguish transient dynamics from steady or asymptotic behavior. Access to larger domains and longer times therefore enables more reliable grain-size distributions, reduces finite-size uncertainty, and strengthens comparisons among theoretical predictions, numerical simulations, and grain-growth behavior observed in real materials [64]. The same acceleration also makes more complex models accessible. As shown by the multicomponent RCD model, each Fourier update becomes more demanding as the number of interacting components increases. By reducing repeated transformed-field work, sp-DBA lowers the computational barrier to studying pattern formation in systems with more interacting species [49] and exploring wider ranges of reaction networks, diffusion parameters, and pattern-forming regimes [65]. In optical communications, simulations of multimode signal transmission likewise become more demanding as the number of modes and operating conditions increases [66]. These examples point to richer models and broader parameter regimes as natural uses of the saved computation. The present results also define the regimes in which sp-DBA is most effective. The largest gains occur when the transformed-field update accounts for a substantial fraction of the computational cost, spectral activity remains concentrated over time, and the active-block ratio remains low. End-to-end acceleration is nevertheless limited by the share of Regular runtime spent in the update stage, as transform operations, communication, and other model steps still contribute to runtime, consistent with Amdahl’s law [67]. A second limitation is that the activation criteria must identify the spectral regions whose updates affect the solution. In the three demonstrations, the dominant changes remain confined to limited spectral regions or distinct bands. More general optical-fiber propagation regimes, however, may involve nonlinear phase evolution with little amplitude change and four-wave mixing that transfers energy across frequencies [68, 69]. Looking ahead, phase-aware criteria, more conservative activation rules, and refresh safeguards offer natural extensions toward regimes with weak amplitude signatures, while dynamic load balancing is the most immediate extension of the distributed implementation. Other directions include block layouts designed for specialized transform hardware and activation rules suited to additional transform families. More broadly, sp-DBA points toward a hardware-aware form of adaptive scientific computing in which high-accuracy numerical methods are paired with 13

controlled adaptive execution to achieve higher performance and broader computational reach.

Acknowledgements We thank Y. Wu for useful discussions.

Author contributions J.J. and Y.X. conceived the study and designed the sp-DBA framework. J.J. implemented the code, performed the simulations, analyzed the data, and drafted the manuscript. P.D. contributed to validation and manuscript review and editing. Y.X. supervised the study and revised the manuscript. All authors reviewed and approved the manuscript.

Competing interests The authors declare no competing interests.

References [1] John P Boyd. Chebyshev and Fourier Spectral Methods. Dover Publications, Mineola, NY, 2nd rev. edition, 2001. [2] Lloyd N. Trefethen. Spectral Methods in MATLAB. Number 10 in Software, Environments, Tools. Society for Industrial and Applied Mathematics, 2000. [3] Claudio Canuto, M. Youssuff Hussaini, Alfio Quarteroni, and Thomas A. Zang. Spectral Methods: Fundamentals in Single Domains. Springer, Berlin, Heidelberg, 2006. [4] James W Cooley and John W Tukey. An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation, 19(90):297–301, 1965. [5] M. Frigo and S.G. Johnson. FFTW: An adaptive software architecture for the FFT. In Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP ’98 (Cat. No.98CH36181), volume 3, pages 1381–1384. IEEE, 1998. [6] NVIDIA. cuFFT Library, 2025. NVIDIA CUDA Documentation; Accessed 17 July 2026. [7] Dmitry Pekurovsky. P3DFFT: A Framework for Parallel Computations of Fourier Transforms in Three Dimensions. SIAM Journal on Scientific Computing, 34(4):C192–C209, 2012. [8] Alan Ayala, Stanimire Tomov, Azzam Haidar, and Jack Dongarra. heFFTe: Highly Efficient FFT for Exascale. In Computational Science—ICCS 2020, volume 12137 of Lecture Notes in Computer Science, pages 262–275, Cham, 2020. Springer International Publishing. [9] NVIDIA. cuFFTMp Documentation, 2025. NVIDIA HPC SDK Documentation; Accessed 17 July 2026. [10] Matti Schneider. A review of nonlinear FFT-based computational homogenization methods. Acta Mechanica, 232(6):2051–2100, 2021. [11] Christian Gierden, Julian Kochmann, Johanna Waimann, Bob Svendsen, and Stefanie Reese. A Review of FE-FFT-Based Two-Scale Methods for Computational Modeling of Microstructure Evolution and Macroscopic Material Behavior. Archives of Computational Methods in Engineering, 29(6):4115–4135, 2022. [12] Claudio Canuto, M. Youssuff Hussaini, Alfio Quarteroni, and Thomas A. Zang. Spectral Methods in Fluid Dynamics. Springer, Berlin, Heidelberg, 1988.

14

[13] G. Kresse and J. Furthmüller. Efficient iterative schemes for Ab Initio total-energy calculations using a plane-wave basis set. Physical Review B, 54(16):11169–11186, 1996. [14] M. D. Feit, J. A. Fleck, and A. Steiger. Solution of the schrödinger equation by a spectral method. Journal of Computational Physics, 47(3):412–433, 1982. [15] Jong Chul Ye. Compressed sensing MRI: A review from signal processing perspective. BMC Biomedical Engineering, 1(1):8, 2019. [16] Ezra Ip and Joseph M. Kahn. Compensation of Dispersion and Nonlinear Impairments Using Digital Backpropagation. Journal of Lightwave Technology, 26(20):3416–3425, 2008. [17] Fatih Yaman and Guifang Li. Nonlinear impairment compensation for polarization-division multiplexed WDM transmission using digital backward propagation. IEEE Photonics Journal, 1(2):144–152, 2009. [18] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations. OpenReview.net, 2021. [19] Wenhui Peng, Shaoxiang Qin, Senwen Yang, Jianchun Wang, Xue Liu, and Liangzhu (Leon) Wang. Fourier neural operator for real-time simulation of 3D dynamic urban microclimate. Building and Environment, 248:111063, 2024. [20] Haitham Hassanieh, Piotr Indyk, Dina Katabi, and Eric Price. Simple and practical algorithm for sparse Fourier transform. In Proceedings of the Forty-Fourth Annual ACM Symposium on Theory of Computing, STOC ’12, pages 1183–1194, New York, NY, USA, 2012. Association for Computing Machinery. [21] Anna C. Gilbert, Piotr Indyk, Mark Iwen, and Ludwig Schmidt. Recent developments in the sparse Fourier transform: A compressed Fourier transform for big data. IEEE Signal Processing Magazine, 31(5):91–100, 2014. [22] Jochen Fröhlich and Kai Schneider. Computation of decaying turbulence in an adaptive wavelet basis. Physica D: Nonlinear Phenomena, 134(3):337–361, 1999. [23] Stephane Mallat. A Wavelet Tour of Signal Processing: The Sparse Way. Academic Press, 3 edition, 2009. [24] Tobias Kaiser, Thorsten Raasch, Joris J.C. Remmers, and Marc G.D. Geers. A waveletenhanced adaptive hierarchical FFT-based approach for the efficient solution of microscale boundary value problems. Computer Methods in Applied Mechanics and Engineering, 409:115959, 2023. [25] Cecilia Pagliantini, Gian Luca Delzanno, and Stefano Markidis. Physics-based adaptivity of a spectral method for the Vlasov–Poisson equations based on the asymmetrically-weighted Hermite expansion in velocity space. Journal of Computational Physics, 488:112252, 2023. [26] Philipp Krah, Thomas Engels, Kai Schneider, and Julius Reiss. Wavelet adaptive proper orthogonal decomposition for large-scale flow data. Advances in Computational Mathematics, 48(2):10, 2022. [27] Alireza Nejadmalayeri, Alexei Vezolainen, Eric Brown-Dymkoski, and Oleg V. Vasilyev. Parallel adaptive wavelet collocation method for PDEs. Journal of Computational Physics, 298:237–253, 2015. 15

[28] Huanxin Lin and Cho-Li Wang. On-GPU thread-data remapping for nested branch divergence. Journal of Parallel and Distributed Computing, 139:75–86, 2020. [29] Albert Segura, Jose Maria Arnau, and Antonio Gonzalez. Irregular accesses reorder unit: Improving GPGPU memory coalescing for graph-based workloads. The Journal of Supercomputing, 79(1):762–787, 2023. [30] Yu-Hsuan Shih, Garrett Wright, Joakim Andén, Johannes P. Blaschke, and Alex H. Barnett. cuFINUFFT: A Load-Balanced GPU Library for General-Purpose Nonuniform FFTs. In 2021 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 688–697. IEEE, 2021. [31] Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale. AMReX: A framework for block-structured adaptive mesh refinement. Journal of Open Source Software, 4(37):1370, 2019. [32] Hsi-Yu Schive, John A. ZuHone, Nathan J. Goldbaum, Matthew J. Turk, Massimo Gaspari, and Chin-Yu Cheng. GAMER-2: A GPU-Accelerated Adaptive Mesh Refinement Code— Accuracy, Performance, and Scalability. Monthly Notices of the Royal Astronomical Society, 481(4):4815–4840, 2018. [33] Andrew Giuliani and Lilia Krivodonova. Adaptive mesh refinement on graphics processing units for applications in gas dynamics. Journal of Computational Physics, 381:67–90, 2019. [34] Ruoyao Zhang and Yang Xia. A dynamic block activation framework for continuum models. Nature Computational Science, 5(4):345–354, 2025. [35] Cody J. Permann, Derek R. Gaston, David Andrš, Robert W. Carlsen, Fande Kong, Alexander D. Lindsay, Jason M. Miller, John W. Peterson, Andrew E. Slaughter, Roy H. Stogner, and Richard C. Martineau. MOOSE: Enabling massively parallel multiphysics simulation. SoftwareX, 11:100430, 2020. [36] K. R. Elder, Mark Katakowski, Mikko Haataja, and Martin Grant. Modeling Elasticity in Crystal Growth. Physical Review Letters, 88(24):245701, 2002. [37] K. R. Elder and Martin Grant. Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals. Physical Review E, 70(5):051605, 2004. [38] Michael Greenwood, Nikolas Provatas, and Jörg Rottler. Free Energy Functionals for Efficient Phase Field Crystal Modeling of Structural Phase Transformations. Physical Review Letters, 105(4):045702, 2010. [39] Michael Greenwood, Jörg Rottler, and Nikolas Provatas. Phase-field-crystal methodology for modeling of structural transformations. Physical Review E, 83(3):031601, 2011. [40] Heike Emmerich, Hartmut Löwen, Raphael Wittkowski, Thomas Gruhn, Gyula I. Tóth, György Tegze, and László Gránásy. Phase-field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: An overview. Advances in Physics, 61(6):665–743, 2012. [41] Yang Xia, Joel M. Berry, and Mikko P. Haataja. Classification and Simulation of Structural Phase Transformation-Induced Interfacial Defects in Group VI Transition-Metal Dichalcogenide Monolayers. Nano Letters, 23(20):9445–9450, 2023.

16

[42] Tatu Pinomaa, Jukka Aho, Jaarli Suviranta, Paul Jreidini, Nikolas Provatas, and Anssi Laukkanen. OpenPFC: An open-source framework for high performance 3D phase field crystal simulations. Modelling and Simulation in Materials Science and Engineering, 32(4):045002, 2024. [43] K. Barmak, E. Eggeling, D. Kinderlehrer, R. Sharp, S. Ta’asan, A.D. Rollett, and K.R. Coffey. Grain growth and the puzzle of its stagnation in thin films: The curious tale of a tail and an ear. Progress in Materials Science, 58(7):987–1055, 2013. [44] R. Backofen, K. Barmak, K.E. Elder, and A. Voigt. Capturing the complex physics behind universal grain size distributions in thin metallic films. Acta Materialia, 64:72–77, 2014. [45] Gabriel Martine La Boissonière and Rustum Choksi. Atom based grain extraction and measurement of geometric properties. Modelling and Simulation in Materials Science and Engineering, 26(3):035001, 2018. [46] Alan M. Turing. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 237(641):37–72, 1952. [47] Nanako Shigesada, Kohkichi Kawasaki, and Ei Teramoto. Spatial segregation of interacting species. Journal of Theoretical Biology, 79(1):83–99, 1979. [48] Dirk Horstmann. Remarks on some Lotka–Volterra type cross-diffusion models. Nonlinear Analysis: Real World Applications, 8(1):90–117, 2007. [49] Edgardo Villar-Sepúlveda, Alan R. Champneys, and Andrew L. Krause. Designing reaction– cross-diffusion systems with Turing and wave instabilities. Journal of Mathematical Biology, 91:37, 2025. [50] Maria Chiara D’Autilia, Ivonne Sgura, and Valeria Simoncini. Matrix-oriented discretization methods for reaction–diffusion PDEs: Comparisons and applications. Computers & Mathematics with Applications, 79(7):2067–2085, 2020. [51] J. Swift and P. C. Hohenberg. Hydrodynamic fluctuations at the convective instability. Physical Review A, 15(1):319–328, 1977. [52] D. J. Richardson, J. M. Fini, and L. E. Nelson. Space-division multiplexing in optical fibres. Nature Photonics, 7(5):354–362, 2013. [53] Antonio Mecozzi, Cristian Antonelli, and Mark Shtaif. Coupled Manakov equations in multimode fibers with strongly coupled groups of modes. Optics Express, 20(21):23436–23441, 2012. [54] Sami Mumtaz, René-Jean Essiambre, and Govind P. Agrawal. Nonlinear propagation in multimode and multicore fibers: Generalization of the Manakov equations. Journal of Lightwave Technology, 31(3):398–406, 2013. [55] Qirui Fan, Gai Zhou, Tao Gui, Chao Lu, and Alan Pak Tao Lau. Advancing theoretical understanding and practical performance of signal processing for nonlinear optical communications through machine learning. Nature Communications, 11(1):3694, 2020. [56] D. Rafique, M. Mussolin, M. Forzati, J. Martensson, M. Chugtai, and A. D. Ellis. Compensation of intra-channel nonlinear fibre impairments using a simplified digital back-propagation algorithm. Optics Express, 19(10):9453–9460, 2011. [57] A. D. Ellis, M. E. McCarthy, M. A. Z. Al Khateeb, M. Sorokina, and N. J. Doran. Performance limits in optical communications due to fiber nonlinearity. Advances in Optics and Photonics, 9(3):429–503, 2017. 17

[58] T. R. Taha and M. J. Ablowitz. Analytical and numerical aspects of certain nonlinear evolution equations. ii. numerical, nonlinear schrödinger equation. Journal of Computational Physics, 55(2):203–230, 1984. [59] Gilbert Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506–517, 1968. [60] O.V. Sinkin, R. Holzlohner, J. Zweck, and C.R. Menyuk. Optimization of the split-step Fourier method in modeling optical-fiber communications systems. Journal of Lightwave Technology, 21(1):61–68, 2003. [61] Moshe Nazarathy and Alex Tolmachev. Digital sub-banding – a signal processing architecture radically improving ofdm coherent optical receivers. In Advanced Photonics Congress, OSA Technical Digest, page SpW2B.1. Optica Publishing Group, 2012. [62] R. Asif, C. Y. Lin, M. Holtmannspoetter, and B. Schmauss. Logarithmic step-size based digital backward propagation in n-channel 112 gbit/s/ch dp-qpsk transmission. In Proceedings of the 13th International Conference on Transparent Optical Networks, page Tu.P.6. IEEE, 2011. [63] Mingtao Xia, Sihong Shao, and Tom Chou. A frequency-dependent p-adaptive technique for spectral methods. Journal of Computational Physics, 446:110627, 2021. [64] Eisuke Miyoshi, Tomohiro Takaki, Munekazu Ohno, Yasushi Shibuta, Shinji Sakane, Takashi Shimokawabe, and Takayuki Aoki. Ultra-large-scale phase-field simulation study of ideal grain growth. npj Computational Materials, 3:25, 2017. [65] Luciano Marcon, Xavier Diego, James Sharpe, and Patrick Müller. High-throughput mathematical analysis identifies Turing networks for patterning with equally diffusing signals. eLife, 5:e14022, 2016. [66] Marius Brehler, Malte Schirwon, Peter M. Krummrich, and Dominik Göddeke. Simulation of nonlinear signal propagation in multimode fibers on multi-GPU systems. Communications in Nonlinear Science and Numerical Simulation, 84:105150, 2020. [67] Gene M. Amdahl. Validity of the single processor approach to achieving large scale computing capabilities. In Proceedings of the April 18–20, 1967, Spring Joint Computer Conference, pages 483–485. ACM Press, 1967. [68] Oleg Sidelnikov, Mikhail Fedoruk, and Stefan Wabnitz. Numerical algorithms for nonlinear propagation in multimode optical fiber communication systems. Optical Fiber Technology, 84:103724, 2024. [69] Peter A. Andrekson. Applications of nonlinear four-wave mixing in optical communication. Optik, 279:170740, 2023. [70] John W. Cahn and John E. Hilliard. Free energy of a nonuniform system. i. interfacial free energy. The Journal of Chemical Physics, 28(2):258–267, 1958. [71] S.M. Cox and P.C. Matthews. Exponential Time Differencing for Stiff Systems. Journal of Computational Physics, 176(2):430–455, 2002. [72] Aly-Khan Kassam and Lloyd N. Trefethen. Fourth-Order Time-Stepping for Stiff PDEs. SIAM Journal on Scientific Computing, 26(4):1214–1233, 2005. [73] Yunfeng Jia, Jianhua Wu, and Hong-Kun Xu. Positive solutions of a Lotka–Volterra competition model with cross-diffusion. Computers & Mathematics with Applications, 68(10):1220–1228, 2014. 18

[74] M. Delgado, M. Montenegro, and A. Suárez. A Lotka–Volterra symbiotic model with cross-diffusion. Journal of Differential Equations, 246(5):2131–2149, 2009. [75] James Ahrens, Berk Geveci, and Charles Law. ParaView: An End-User Tool for Large-Data Visualization. In Charles D. Hansen and Chris R. Johnson, editors, Visualization Handbook, pages 717–731. Elsevier, 2005. [76] Alexander Stukowski. Visualization and analysis of atomistic simulation data with OVITO— the open visualization tool. Modelling and Simulation in Materials Science and Engineering, 18(1):015012, 2010.

19

Supplementary Information sp-DBA framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 S1.1 Boundary-aware transition smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 S1.2 Neighbor activation in half-spectrum FFT layouts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 S1.3 Dynamic threshold control and periodic refresh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 S2 sp-DBA integration with transform-domain solvers. . . . . . . . . . . . . . . . . . . . . . . . . . . .25 S2.1 Integration procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 S2.2 Implementation overhead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 S2.3 Comparison with adaptive spectral and wavelet methods . . . . . . . . . . . . . . . . . . . . . . . . . 27 S3 XPFC polycrystalline coarsening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 S3.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 S3.2 Numerical scheme and simulation configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 S3.3 Grain reconstruction and system-size comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 S4 Reaction–cross-diffusion pattern formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 S4.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 S4.2 Numerical scheme and simulation configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 S4.3 Matrix eigendecomposition for GPU acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 S5 Digital backpropagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 S5.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 S5.2 Split-step Fourier implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 S5.3 Simulation configuration and recovery metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 S6 Numerical validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 S6.1 Conservation and free-energy evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 S6.2 Spectral accuracy and numerical error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 S6.3 Dealiasing compatibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 S7 Single-GPU performance evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 S7.1 PFC cost decomposition and acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 S7.2 DBP cost decomposition and acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 S8 Multi-GPU sp-DBA with distributed FFTs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .42 S8.1 Distributed FFT workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 S8.2 Strong and weak scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 S8.3 CUDA-aware MPI communication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 S8.4 Activation stability under distributed execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 S9 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 S1

20

S1 sp-DBA framework sp-DBA uses boundary-aware transition smoothing, half-spectrum neighbor activation, dynamic threshold control, and periodic refresh to manage block states and update boundaries in transformed fields. Hereafter, Regular refers to the full-update mode with sp-DBA disabled.

S1.1 Boundary-aware transition smoothing When the solver update is skipped in frozen blocks, their transformed-field values are retained from the previous step. A moving boundary can therefore appear between updated and skipped regions in the one-step update increment. To reduce this boundary discontinuity, sp-DBA applies a smooth weight to the newly applied increment, q̂n+1 (k) = q̂n (k) + σ(k)∆q̂raw (k),

(S.1)

where q̂ is a generic transformed variable, ∆q̂raw is the raw update increment, and σ(k) is the transition weight. The weight is applied only to the increment, so the stored values in frozen blocks are not damped by the smoothing operation [1, 3]. For a block with edge length B, the one-dimensional taper is defined as πi 1 1 − cos 2 B−1 

w(i) =





,

i = 0, 1, . . . , B − 1.

(S.2)

The current block state and neighboring block states are loaded into shared memory, and the local transition weight is assembled from the frozen-neighbor pattern. In two dimensions, let χL , χR , χD , χU ∈ {0, 1} indicate whether the left, right, lower and upper neighboring blocks are frozen. The local weight is then σ2D (i, j; t) = [w(i)]χL [w(B − 1 − i)]χR [w(j)]χD [w(B − 1 − j)]χU .

(S.3)

The same construction extends to three dimensions by adding the front and back factors. Only the one-dimensional taper and block-state metadata are stored, while the multidimensional weight is assembled locally during the update.

21

Boundary-aware transition update Active Transition

1 Roll-off factor

Frozen 0

State-aware by shared memory

State-aware access Figure S1: Boundary-aware transition smoothing. Neighboring block states are loaded through shared memory to form a local roll-off weight. Active blocks receive the full update increment, transition blocks receive a weighted increment, and frozen blocks retain their previous transformed-field values.

no smooth

smooth

b

sp-DBA state

spectrum

a

Figure S2: Effect of transition smoothing. a, Fourier-amplitude maps and sp-DBA block states without and with smoothing. b, Total free-energy evolution for Regular, sp-DBA without smoothing, and sp-DBA with smoothing.

Smoothing leaves the block pattern largely unchanged and brings the free-energy trajectory closer to that of Regular (Fig. S2).

22

S1.2 Neighbor activation in half-spectrum FFT layouts For real-valued fields, the PFC and RCD solvers store only the nonredundant half of Fourier space in the real-to-complex FFT layout. sp-DBA applies neighbor activation on this stored half-spectrum while preserving the conjugate symmetry required by real-valued fields. If a requested neighbor lies inside the stored half-spectrum, it is activated directly. If the requested neighbor lies in the omitted half, sp-DBA activates its conjugate counterpart using q̂(−k) = q̂ ∗ (k).

a

b

transition block

active block

(S.4)

DC -5

-4

-3

-2

-1

0

1

2

3

4

5

6

active block

4

5

-6

3

mark neighboring blocks

yes

no

wake directly

conjugate remap

-4

-3

-2

-1

0

1

2

Periodic boundary

stored half?

-5

wake symmetric partner

conjugate counterpart (no stored)

half spectrum (actual stored)

Figure S3: Half-spectrum neighbor activation. a, Neighbor activation on the stored half-spectrum. Neighbor blocks outside the stored half are activated through their conjugate counterparts. b, Simplified decision rule. Neighbors inside the stored half are activated directly; neighbors outside it are mapped by conjugate symmetry, and the symmetric partner on the kx = 0 column is activated when required.

The kx = 0 column is treated as a symmetry boundary. When conjugate remapping reaches this column, the corresponding symmetric partner is also activated. This remapping affects only block-state propagation and does not modify the stored Fourier coefficients.

S1.3 Dynamic threshold control and periodic refresh Dynamic threshold. Dynamic threshold control adjusts the activation threshold to maintain a target average active-block ratio. At each control interval, the measured average active-block ratio r, including active and transition blocks, is compared with the target ratio rtar . The threshold is increased when r > rtar and decreased when r < rtar . For PFC, deviations outside a relative deadband η use factors 1 + η or 1 − η, with factor-of-two corrections when r/rtar > 2 or r/rtar < 0.5.

23

a

40962

b

10242

c

20482

d

40962

e

81922

Figure S4: Dynamic threshold control in the XPFC application. a, Average active-block ratio on a 40962 grid. The upper plot shows runs with target ratios of 5%, 10%, 15%, and 20%; the lower plot shows fixed-threshold behavior during the first 10,000 steps, which are excluded from the reported performance statistics. b–e, Runs with dynamic threshold control on 10242 , 20482 , 40962 , and 81922 grids. After the initial full-activation stage, the measured active-block ratio stays close to its target.

In the reported PFC runs, control begins after the initial 10,000-step full-activation stage, and each ratio check is followed by one full-activation step. Periodic refresh. Periodic refresh activates all blocks every Frefresh steps to limit long-time drift. Between refresh events, blocks selected by the activation rule are updated at the base time step; for example, ∆t = 0.01 and Frefresh = 50 give a refresh interval of 0.5 time units.

24

a

∆t Active Block

Frozen Block

t Hold last (skipped)

Update

b

Hold last (skipped)

Full activation Frefresh·∆t

t Full activation

c

Figure S5: Periodic refresh. a, Active blocks are updated at each step ∆t. Frozen blocks retain their previous state between local reactivation events or full-activation events separated by Frefresh ∆t. b, Total free-energy evolution for Regular and sp-DBA with different refresh intervals. The inset shows the mean density ψ̄. c, Average active-block ratio for the same refresh intervals.

Shorter refresh intervals bring the free-energy curves closer to the Regular trajectory, while the average active-block ratio returns to its low baseline after each refresh (Fig. S5). The mean density is unchanged over the tested intervals.

S2 sp-DBA integration with transform-domain solvers Figure S6 summarizes the shared sp-DBA interface used in the three solvers.

25

Figure S6: Integration of sp-DBA into the three solver workflows. DBP connects sp-DBA to the two frequency-domain linear half-steps of the one-dimensional full-spectrum split-step solver, while PFC and RCD connect it to the transformed-field update of two- and three-dimensional pseudo-spectral solvers. Regular and sp-DBA call the same model-specific update_core; sp-DBA controls which blocks execute it and how transition increments are weighted. The shared interface initializes block states, dispatches the update, applies transition handling, updates neighboring states, performs periodic refresh, and releases device arrays.

S2.1 Integration procedure Each solver exposes its transformed-field update as update_core. Regular calls this routine over the full transformed field, whereas sp-DBA calls it only for active and transition blocks. Frozen blocks return early, and transition blocks use the weighted increment defined in Supplementary Section S1. The solver invokes spdba_apply(...) at this stage. With mode=0, it dispatches the full update; with mode=1, it checks block states, calls update_core where required, updates neighboring states, and applies threshold or refresh control. Each model adapter supplies the layout, indexing, launch configuration, and model parameters required by the shared interface. It binds the model-specific update_core to the generic kernel for the one-dimensional full-spectrum DBP layout and the two- and three-dimensional half-spectrum PFC and RCD layouts. Pseudocode for integration of sp-DBA into a solver readInput(...) allocate model fields construct spDBA parameters, state arrays, and model adapter if (params.mode == 1) spdba_init(...) for each simulation step model_pre_update(...) construct spDBA_step(...)

// model terms, transforms, preprocessing

26

spdba_apply(...) model_post_update(...)

// mode=0 full update; mode=1 sp-DBA // inverse transforms, normalization, output

if (params.mode == 1) spdba_destroy(...)

The generic kernels construct block-local indices and call the model-provided update_core for the SPDBA_1D_FULL, SPDBA_2D_HALF, and SPDBA_3D_HALF layouts.

S2.2 Implementation overhead Implementation overhead is measured in nonblank lines of code (LOC), following Ref. [34], with separate counts for the shared core, model adapter, and solver integration. The shared spDBA/ core contains 696 LOC and implements block partitioning, state management, transition smoothing, neighbor activation, threshold control, and refresh; it is excluded from the counts and percentages for each solver. Simple changes to each solver include exposing update_core, binding the sp-DBA state, reading parameters, connecting diagnostic output, and calling the shared interface. Input/output code, scripts, documentation, tests, generated files, and historical files are excluded. Table S1: Implementation overhead for integrating sp-DBA into the three solvers. Application Regular core Solver LOC Solver % Adapter Combined LOC Combined % API calls PFC RCD DBP

2162 1284 1513

178 211 159

8.23% 16.43% 10.51%

117 145 97

295 356 256

13.64% 27.73% 16.92%

9 11 19

The simple changes add 159–211 nonblank LOC, corresponding to about 8–16% of each Regular solver core. Including the model adapters increases the total to 256–356 LOC, or about 13–28% of the original core sizes. RCD has the largest percentage because its integration includes a three-dimensional half-spectrum layout, multicomponent fields, eigendecomposition, legacy update paths, and shared-memory launch configuration.

S2.3 Comparison with adaptive spectral and wavelet methods Table S2 compares how representative adaptive spectral, wavelet, and sp-DBA workflows introduce adaptivity, focusing on changes to the numerical representation, data structures, and solver implementation. Runtime and accuracy are not compared across methods. Representative adaptive spectral methods modify the approximation space during the calculation. Frequency-dependent p-adaptivity adjusts the expansion order and reconstructs the solution on a new set of collocation points, while adaptive Hermite methods update basis parameters and map the numerical state between successive approximation spaces [63, 25]. Adaptive wavelet workflows introduce a hierarchy of resolution levels. Wavelet-enhanced FFT methods refine the computational grid and reformulate the Fourier-space operator for the wavelet discretization, while parallel implementations use hierarchical storage, refinement, coarsening, reconstruction or ghost regions, and, in distributed calculations, block migration and dynamic load balancing [24, 27, 26]. sp-DBA instead retains the transformed representation, regular array layout, transform operations, and model-specific update core, while adding block state arrays and activation control around the existing update. Table S2 summarizes the contrast.

27

Table S2: Representative implementation paths for adaptive spectral and wavelet workflows and sp-DBA. Adaptive strategy

Quantity changed during the calculation

Additional implementation path

Adaptive spectral basis

Expansion order or basis parameters Active resolution levels, grid points or wavelet coefficients

Reconstruct or project the numerical state between successive approximation spaces Maintain hierarchical grid data and manage refinement, coarsening, reconstruction or ghost regions, and, where required, migration and load balancing Maintain block states and activation control around the retained update core; the representation, regular data layout, and transform operations remain in place

Wavelet and multiresolution workflow sp-DBA

Blocks that execute the existing transformed-field update

Implementation requirements vary among adaptive spectral and wavelet methods. Table S1 reports the added code for the three sp-DBA integrations and is not a cross-method LOC comparison.

S3 XPFC polycrystalline coarsening The XPFC tests use a scalar pseudo-spectral model for polycrystalline coarsening.

S3.1 Governing equations We use a dimensionless XPFC formulation [38, 41]. The order parameter ψ(r) describes deviations from a reference density and evolves according to a free energy with local polynomial terms and a nonlocal two-point correlation term, Z

F [ψ] =

1 ω u 1 dr ψ 2 − ψ 3 + ψ 4 − 2 6 12 2 



ZZ

dr dr ′ ψ(r)C2 (|r − r ′ |)ψ(r ′ ).

(S.5)

The local terms define the bulk free energy and stabilize the density amplitude, while the nonlocal term favors density modulations selected by the correlation kernel C2 . The conserved dynamics follow the Cahn–Hilliard form [70], ∂ψ δF = ∇2 , (S.6) ∂t δψ which gives ∂ψ = ∇2 ψ − ∂t 

Z

dr ′ C2 (|r − r ′ |)ψ(r ′ ) −

ω 2 u 3 ψ + ψ . 2 3 

(S.7)

The outer Laplacian enforces conservation of the mean density. In Fourier space, the convolution with C2 becomes multiplication by Ĉ2 (k). For the twodimensional hexagonal XPFC setting used here, the reciprocal-space kernel is represented by a single Gaussian peak, " # " # (k − q0 )2 σ 2 q02 Ĉ2 (k) = exp − exp − , (S.8) 2α2 2ρβ with ρ = β = 1 in the production coarsening runs. The first exponential defines the structural peak around the preferred wavenumber q0 , and the second exponential acts as a Debye–Waller attenuation factor. The pseudo-spectral implementation evaluates the local nonlinear terms in real space and the correlation term in Fourier space. In polycrystalline states, different grain orientations broaden the Bragg peaks into ring-like Fourier-amplitude regions, which are used to assign the PFC block states.

28

S3.2 Numerical scheme and simulation configuration Time-stepping method. The simulations are advanced in Fourier space with the fourth-order multistep exponential time-differencing scheme (ETD4) of Cox and Matthews [71]. ETD4 advances the linear term through the exponential propagator eL∆t and evaluates the nonlinear contribution using its current value and three previous Fourier-space values. Writing the Fourier-space evolution as dψb \ = L(k)ψb + N (ψ), dt

(S.9)

the update is 







b n − 2N b n−1 + N b n−2 Φ3 b n Φ1 + N bn−N b n−1 Φ2 + N ψb n+1 = ez ψb n + N 



b n − 3N b n−1 + 3N b n−2 − N b n−3 Φ4 , + N

(S.10)

bn =N \ where N (ψ n ), z = L(k)∆t, and

Φ1 = ∆t ϕ1 (z), Φ2 = ∆t ϕ2 (z),   1 Φ3 = ∆t ϕ2 (z) + ϕ3 (z) , 2   1 Φ4 = ∆t ϕ2 (z) + ϕ3 (z) + ϕ4 (z) . 3

(S.11)

The coefficient functions are defined by ϕℓ (z) =

∞ X

zj , (ℓ + j)! j=0

ℓ = 1, . . . , 4.

(S.12)

For small z, the coefficient functions are evaluated using their Taylor expansions to avoid cancellation and apparent division by zero [71, 72]. In the conserved PFC implementation, \ (ψ), and the corresponding factor −k 2 ∆t is the outer Fourier-space Laplacian is included in N incorporated into the precomputed nonlinear-update coefficients. At initialization, the three nonlinear-history arrays are set to the initial nonlinear term and shifted after each time step. Simulation and initialization parameters. The XPFC simulations are performed on a two-dimensional periodic square domain. Performance tests use grid sizes from 5122 to 163842 . The layered main-text visualization uses a 20482 production run, whereas the main-text GSD and free-energy comparisons use 81922 production runs. The representative production configuration uses ∆x = ∆y = 0.7854, ∆t = 0.01, and 1, 000, 000 time steps, corresponding to a final dimensionless time of t = 10, 000. The XPFC evolution parameters are q0 = 1.0, σ = 0.02, α = 0.1, ω = 1.0, and u = 1.0. Here, q0 , σ, and α parameterize the reciprocal-space correlation kernel, while ω and u define the local nonlinear terms in the free energy. The polycrystalline initial condition is constructed around a mean density ψc = 0.12 by adding randomly oriented crystalline seeds with density modulation amplitude Ah = 0.25. The seed centers and orientations are generated from fixed random seeds and are smoothed near seed boundaries before the initial growth stage. For the sp-DBA runs, the Fourier block size is B = 16. The initial threshold is 5 × 10−5 , boundary-aware transition smoothing and dynamic threshold control are enabled, the adjustment factor is 0.05, and the control interval is 50 steps. The target average active-block ratio ranges from 5% to 20%. During the first 10,000 steps, all Fourier blocks are kept active to avoid premature freezing during the rapidly evolving initial transient. The late-stage analysis begins after this full-activation period. 29

Fixed-threshold sensitivity is also tested on a 40962 grid using the same simulation and sp-DBA configuration. Only the threshold policy changes, with fixed values ϵT = 1 × 10−4 , 3 × 10−5 , 4 × 10−5 , and 5 × 10−5 . Across these values, the normalized grain-size distribution and grain-growth trend agree with Regular (Fig. S7). The main text uses dynamic threshold control because it keeps the average active-block ratio more stable across coarsening stages and grid sizes.

a

b

Figure S7: Fixed-threshold sensitivity in the XPFC application. a, Mean squared equivalent grain 2 diameter ⟨Deq ⟩ as a function of sampled coarsening time for Regular and sp-DBA with fixed threshold values ϵT on a 40962 grid. b, Normalized grain-size distribution p(x) with x = Deq /⟨Deq ⟩ for the same fixed threshold settings. The gray band indicates the Regular result range.

Layered correspondence. Regular t = 10

sp-DBA t = 100

sp-DBA t = 10,000

Real space ψ

Fourier space Amplitude

sp-DBA state

Figure S8: Layered XPFC coarsening under sp-DBA on a 2048 × 2048 grid. The top, middle, and bottom rows show the real-space density field ψ, Fourier amplitude, and sp-DBA block state, respectively. The average active-block ratio is about 15%.

30

The real-space density evolves through grain-boundary migration and defect elimination, while the Fourier amplitude and active blocks are concentrated near the dominant Bragg ring.

S3.3 Grain reconstruction and system-size comparison Grain statistics are reconstructed from detected density peaks following Ref. [45]. Because the labels contain post-processing noise, the late-stage scaling plots are treated as qualitative agreement checks [44]. Extraction procedure. The continuous density field is first converted into an atom-like point set by detecting local density maxima under periodic boundary conditions. For each detected peak, a first-neighbor shell is constructed within a cutoff radius, and a local bond-orientational descriptor is evaluated. For the two-dimensional hexagonal lattice, the orientation is represented by the phase of the sixfold bond-order parameter. Atoms with sufficiently ordered local environments are selected as grain cores, while atoms with weak orientational order or strong distortion are treated as boundary or defect candidates. Neighboring core atoms are merged into the same grain when their misorientation remains below the prescribed large-angle boundary threshold. Excluded boundary atoms are then reassigned by neighborhood voting and proximity checks. Very small fragments, enclosed islands and low-misorientation artifacts are removed or merged. The cleaned labels and a local atomic-area approximation are then used to evaluate grain area, equivalent diameter, GSD and coarsening statistics in a form comparable with standard grain-growth analysis [44]. Extraction settings and cross-size comparison. The same reconstruction workflow is applied to the sampled OVITO dump frames from the Regular run and from sp-DBA runs with average active-block ratios of 5%, 10%, 15%, and 20%. The fixed reconstruction parameters are rcut = 1.18, θcut = 3.5◦ , ψkeep = 0.55, and a minimum grain core size of 80 atoms. Boundary atoms are reassigned using six voting passes with an assignment radius factor of 1.3, and small-fragment cleanup uses a merge factor of 3.0. For the normalized GSD, grain diameters are rescaled frame by frame as x = Deq /⟨Deq ⟩, and the steady distribution is constructed from the last 50% of sampled frames. The histogram uses 22, 28, 34, and 40 bins for 10242 , 20482 , 40962 , and 81922 , respectively, with the plotted range limited to 0 < x < 3.5. The supplementary scaling plots use sampled coarsening time t and 2 ⟩ ∼ t0.3 for visual comparison. show the common guide line ⟨Deq

31

a

b

10242

e

c

20482

f

40962

g

d

81922

h

Figure S9: Supplementary GSD and late-stage scaling comparisons across XPFC sizes. a–d, Steady normalized GSD comparisons for the 10242 , 20482 , 40962 , and 81922 production runs. e–h, Corresponding 2 sampled late-stage scaling plots for the same four sizes, shown with the common guide line ⟨Deq ⟩ ∼ t0.3 .

The GSD collapse becomes clearer with increasing system size, with the largest systems giving the smoothest late-stage distributions (Fig. S9).

S4 Reaction–cross-diffusion pattern formation The RCD tests use a constant-coefficient multicomponent model.

S4.1 Governing equations The general multicomponent reaction–cross-diffusion system is written as [48, 73] ∂u = ∇ · (D(u)∇u) + R(u), ∂t

(S.13)

where u = (u1 , . . . , uM )T is the concentration vector, D(u) is the diffusion matrix, and R(u) is the reaction term. In the reported tests, the model is reduced to the constant-coefficient form ∂u = D∇2 u + R(u). ∂t

(S.14)

The diffusion matrix is chosen as a symmetric circulant cross-diffusion matrix over the component index, 

Ddiag + 2Dcross −Dcross  −Dcross Ddiag + 2Dcross 

  D=   

0 −Dcross

0 .. .

−Dcross .. .

Ddiag + 2Dcross .. .

−Dcross

0

···

··· ···

−Dcross 0 .. .

··· .. . −Dcross −Dcross Ddiag + 2Dcross

      . (S.15)   

For M = 2, the two periodic neighbors coincide, and the single off-diagonal entry is −Dcross . Each component is coupled to its two neighboring components under periodic indexing. Because D is constant, the optional eigendecomposition described in Supplementary Section S4.3 is computed once before time stepping [50].

32

The nonlinear reaction term is taken in generalized Lotka–Volterra form [48, 73, 74], 



Ri (u) = ui ri − ai ui −

X

αij uj  .

(S.16)

j̸=i

Here, ri is the intrinsic growth rate, ai is the self-limitation coefficient, and αij describes the interaction from component j to component i. In the reported tests, ri = 1, ai = 1, αi,i+1 = αi,i−1 = β + α, and αij = β for all other off-diagonal entries, again with periodic indexing over components.

S4.2 Numerical scheme and simulation configuration Time stepping method. Time integration uses a semi-implicit Fourier update. The linear cross-diffusion term is treated implicitly, and the nonlinear reaction term is evaluated explicitly. For each stored Fourier mode k, 

û n+1 (k) = I + ∆t k 2 D

−1 h

i

n b û n (k) + ∆t R(u )(k) .

(S.17)

Thus, in Regular, each stored Fourier mode requires an M × M matrix operation. Simulation parameters. The tested configurations use periodic three-dimensional domains with spatial resolutions of 1283 and 2563 , while the number of components is varied from M = 2 to M = 64. The spatial spacing is ∆x = 0.5, the time step is ∆t = 0.01, and the total number of steps is 20,000. The constant diffusion matrix uses Ddiag = 0.5 and Dcross = 0.05. The generalized Lotka–Volterra reaction uses α = 1.0 for the enhanced neighboring component interaction and β = 0.3 for the remaining intercomponent competition. The initial condition is a nearly homogeneous state with random perturbation amplitude 0.1. For the sp-DBA runs, the block size is B = 4 for 1283 and B = 8 for 2563 , the initial threshold is 10−9 , and dynamic threshold control is enabled. The adjustment factor is 0.1, the control interval is 50 steps, and the target average active-block ratio ranges from 3% to 20%. For RCD block selection, the activation criterion is computed from all components rather than from the single component shown in visualizations. At each stored Fourier mode k, the implementation uses the component RMS Fourier amplitude A(k) =

M 1 X |ûm (k)|2 M m=1

!1/2

.

(S.18)

A Fourier block Bj is activated when at least one stored mode satisfies A(k) > εRCD , and the resulting state is shared by all M components in that block. The eigendecomposition of D is used only within the update core and does not enter the activation criterion. Layered correspondence.

33

step = 1,000

step = 4,000

step = 20,000

Real space component 1

Fourier space Amplitude

sp-DBA state

Figure S10: Layered RCD pattern formation under sp-DBA on a 256 × 256 × 256 grid with M = 8. The top, middle, and bottom rows show one real-space component, Fourier amplitude, and sp-DBA block state, respectively. The average active-block ratio is about 3%.

The Fourier amplitude is localized in a compact low-wavenumber region, and the sp-DBA active blocks follow this region (Fig. S10).

S4.3 Matrix eigendecomposition for GPU acceleration For the constant-coefficient tests, the diffusion matrix is diagonalized once before time stepping, D = VΛV−1 . Then



(S.19) 

I + ∆t k 2 D = V I + ∆t k 2 Λ V−1 ,

(S.20)

and the update becomes 

û n+1 (k) = V I + ∆t k 2 Λ

−1

h

i

n b V−1 û n (k) + ∆t R(u )(k) .

(S.21)

This form replaces the dense matrix inverse with two constant basis transforms and one diagonal scaling that depends on k.

34

block

Register

Shared Memory

Global Memory (HBM)

Figure S11: Hardware-oriented diagonalized matrix update in the RCD solver. The gray path shows the direct matrix update, and the red path shows the reordered diagonalized implementation. The diagonal factor is held in registers, while the constant basis matrices are cached in shared memory for reuse within a block.

Eigendecomposition reduces the cost of the coupled update at each stored Fourier mode, whereas sp-DBA reduces the number of Fourier blocks where that update is executed. The two optimizations therefore act on different factors of the RCD update cost. For the timing comparisons in the main text, Regular and sp-DBA use the direct coupled matrix update, whereas Regular + Eigen and sp-DBA + Eigen use the diagonalized update. The 28.1× value at M = 48 in Fig. 3d compares sp-DBA + Eigen with Regular using the direct matrix update; Fig. 5d reports all four configurations to separate the two contributions.

S5 Digital backpropagation The DBP tests use a controlled split-step signal-recovery calculation motivated by space-division multiplexed optical communication [16, 17].

S5.1 Governing equations Signal propagation is described by a simplified multimode model with linear mode coupling [53, 54]. For the complex field envelope vector E = [E1 , E2 , . . . , EM ]T , the propagation equation is ∂E = (D̂ + N̂ )E. ∂z

(S.22)

The linear and nonlinear operators are α β2 ∂ 2 D̂ = − I − i I + iK, 2 2 ∂t2

8 N̂ = iγ 9

M X

! 2

|Em |

I.

(S.23)

m=1

Here, α is the attenuation coefficient, β2 is the group-velocity dispersion coefficient, γ is the effective Kerr coefficient, and K is a constant Hermitian mode-coupling matrix. The attenuation and dispersion terms act independently on each mode, while K mixes the modal components through fixed linear coupling. The nonlinear operator gives the Manakov-averaged Kerr phase rotation driven by the total modal power and is evaluated pointwise in the time domain. In the frequency domain, the dispersion part of D̂ is diagonal over frequency bins, while mode coupling acts as a constant matrix on the modal vector at each frequency. DBP approximates inverse propagation by reversing the dispersion, attenuation, and nonlinear parameters and by applying the inverse mode-coupling propagator [16, 17].

35

S5.2 Split-step Fourier implementation Propagation is evaluated with a symmetric split-step Fourier method (SSFM) [59, 60, 68]. Over a step dz, Z z+dz

dz E(z + dz, t) ≈ exp D̂ exp 2 



! ′

N̂ (z ) dz

′

z

dz exp D̂ E(z, t). 2 



(S.24)

Each propagation step contains two frequency-domain linear half-steps and one time-domain nonlinear step. With both dispersion and mode coupling included in the linear operator, each linear half-step is factorized as 

exp

dz dz exp D̂ ≈ exp iK 2 4 







α β2 − I + i ω2I 2 2



dz dz exp iK , 2 4 





(S.25)

where ω is the angular frequency of the Fourier bin. The quarter-step coupling propagator is Uc (δ) = exp(iδK),

δ=

dz . 4

(S.26)

Because K is Hermitian and constant, it is diagonalized once as K = V ΛV † ,

Uc (δ) = V exp(iδΛ)V † .

(S.27)

Forward and backward coupling actions are obtained by using opposite signs of δ, giving consistent discrete coupling propagators. For the coupled-mode tests, K is generated once using std::mt19937 with seed 1234. Its diagonal entries are real samples from N (0, 1); the real and imaginary parts of each uppertriangular off-diagonal entry are sampled independently from the same distribution, and the lower triangle is set by Hermitian symmetry. The matrix is normalized such that M 1 X |Kij |2 = 1. M 2 i,j=1

(S.28)

sp-DBA is applied only to these linear half-steps. The FFT and IFFT calls, the time-domain nonlinear step, and the split-step ordering are unchanged.

S5.3 Simulation configuration and recovery metrics Simulation parameters. All propagation and signal parameters are expressed in the simulation units used throughout the DBP calculations. The FFT size is varied from Nt = 8192 to 262144, and the number of coupled components is varied from M = 2 to M = 64 in the scaling tests. The propagation configuration uses a fiber length of 500 km, a step size of dz = 0.1 km, and an output interval of 200 steps. The propagation parameters are β2 = −21, α = 0.1 in the exponential attenuation convention used by the simulation, and γ = 1.3. For the sp-DBA runs, the frequency block size is B = 32, the threshold is initialized at 10−4 , and dynamic threshold control is enabled. Unless otherwise stated, the default performance setting targets a 10% average active-block ratio. Additional accuracy and sensitivity tests use target average active-block ratios of 5%, 10%, 15%, and 20%. The threshold control interval is 50 steps, and the dynamic adjustment factor is 0.05. Reported active-block ratios are averages over the recovery run. Deterministic attenuation changes the overall frequency-domain amplitude scale, so the dynamic threshold is rescaled by the running signal level before activation. Initial signal configuration. 36

The transmitted field is initialized as a sum of three Gaussian wave packets with different carrier frequency offsets, 3 X

#

"

(t − tq )2 exp(iΩq t) . Em (0, t) = A0 exp − 2w2 q=1

(S.29)

The same initial waveform is assigned to all modal components. The three packets use temporal offsets tq = (−20, 0, 20) and carrier offsets Ωq = 2π

Nt (−0.15, 0, 0.15), T

(S.30) p

where T = 1000 is the temporal window. Each packet has amplitude A0 = P0 /3, where P0 = 1000 is the prescribed total peak-power parameter, and w = 10 is the prescribed pulse width. Frequency-block activation. For block selection, sp-DBA uses the total modal spectral power P (ω) =

M X

|Ẽm (ω)|2 .

(S.31)

m=1

A frequency block is activated when its spectral power exceeds the current threshold after the attenuation rescaling described above, and the resulting state is applied to the full M -component frequency-domain vector in that block. Reconstruction metrics. Let Etx (t) denote the transmitted field and Ep (t) the recovered field for p ∈ {reg, spDBA}. For each recovered field, the normalized mean-squared error is evaluated against the transmitted reference field, 2 t ∥Ep (t) − Etx (t)∥2 P 2 t ∥Etx (t)∥2

P

NMSEp (dB) = 10 log10

!

.

(S.32)

Here, the sum is taken over the sampled time points and ∥ · ∥2 denotes the Euclidean norm over the modal components. The field correlation coefficient is P

Corrp = 

† t Ep (t) Etx (t)

2 t ∥Ep (t)∥2

P

1/2 P

2 t ∥Etx (t)∥2

1/2 ,

(S.33)

where (·)† denotes the Hermitian transpose.

a

Spectrum evolution Frequency

b

Recovery illustration

c

Activation evolution Frequency

Amplitude

sp-DBA state

Recovery distance

Recovery distance

Figure S12: Spectral evolution and sp-DBA activation in DBP. a, Spectrum versus recovery distance and frequency; color shows spectral amplitude. b, Schematic signal recovery and corresponding sp-DBA block states. c, Block-state evolution along the recovery distance. Low, intermediate, and high values correspond to frozen, transition, and active regions, respectively.

37

The dominant spectral bands remain localized over much of the recovery distance, and the sp-DBA activation pattern follows their evolution (Fig. S12).

S6 Numerical validation S6.1 Conservation and free-energy evolution For the XPFC application, we use the conserved mean-density mode and the free-energy trend to check whether block activation introduces systematic drift. The block containing the zerofrequency mode is kept active throughout the run, so the conserved mode is advanced as in Regular.

a

10242

b

20482

c

40962

d

81922

Figure S13: Free-energy and mean-density agreement in the XPFC application under sp-DBA. a–d, Total free-energy evolution for Regular and sp-DBA on 10242 , 20482 , 40962 , and 81922 grids. The insets show the conserved mean density ψ̄ for the corresponding runs.

Across the tested grid sizes, the mean density agrees with Regular, with no systematic increase in the free-energy discrepancy (Fig. S13).

38

S6.2 Spectral accuracy and numerical error

a

b

Frozen region: values fall below 10-16

Figure S14: Spectral agreement and numerical error under sp-DBA. a, Radial structure factor corresponding to the main-text comparison. Away from the principal peak region, the spectrum decays rapidly and the tail reaches the machine-precision level. b, Relative error after enabling sp-DBA.

Outside the principal peak region, the radial structure factor decays to approximately 10−16 , while the relative discrepancy after enabling sp-DBA remains about 10−7 (Fig. S14). This residual is consistent with the truncation-type error expected when spectrally small regions are retained from previous steps [1].

39

S6.3 Dealiasing compatibility sp-DBA

sp-DBA + dealiasing

Fourier space

Dealiasing areas

sp-DBA state

Figure S15: Compatibility of sp-DBA with dealiasing. The two columns compare sp-DBA without and with dealiasing. In each column, the upper row shows the corresponding spectrum and the lower row shows the sp-DBA block activation state.

Pseudo-spectral evaluation of nonlinear terms can produce aliasing when unresolved highfrequency interactions fold back into the represented spectral band [1]. Dealiasing first defines the retained spectral region, within which sp-DBA applies block activation (Fig. S15).

S7 Single-GPU performance evaluation All single-GPU performance measurements were obtained on one NVIDIA A800-SXM4 GPU using the same CUDA version, warm-up procedure, repeat count and GPU-synchronized wallclock timing protocol listed for the distributed tests in Table S3, with MPI disabled. Reported timings exclude file I/O and post-processing. For runs in which communication is absent or not dominant, the measured end-to-end acceleration is approximated by Stotal =

Treg Tsp-DBA

≈

Treg , TFFT + TNL + αTupdate + Toverhead

(S.34)

Here, TFFT , TNL , and Tupdate denote the transform, retained model, and Regular update costs, respectively; α is the average active-block ratio, and Toverhead includes state evaluation, transition handling, refresh, and control operations. The unchanged workflow components impose the corresponding Amdahl-type bound [67]. A simpler estimate based on the active-block ratio is Sestimate ≈

TFFT + TNL + Tupdate , TFFT + TNL + αTupdate 40

(S.35)

where sp-DBA control overhead is omitted. Departures from this estimate arise from the omitted control overhead and hardware-dependent workload, memory, and execution effects. For the RCD comparison in Fig. 5e of the main text, Eq. (S.35) serves as a block-activation estimate, while the measured curve also includes eigendecomposition.

S7.1 PFC cost decomposition and acceleration

a

e

Ratioavg= 5%

b

f

Ratioavg= 10%

c

g

Ratioavg= 15%

d

Ratioavg= 20%

h

Figure S16: Additional end-to-end performance analysis for the XPFC application. a–d, Total wall-clock cost decomposition for Regular and sp-DBA at average active-block ratios of 5%, 10%, 15%, and 20%, respectively. e–h, Comparison between the estimate from Eq. (S.35) and the measured end-to-end acceleration for the same four active-block ratios.

Most of the XPFC reduction comes from the transformed-field update, while the transform and retained model costs remain nearly unchanged (Fig. S16a–d). At larger grid sizes, the increased update workload better amortizes the block-control overhead, bringing the measured acceleration closer to the active-block-ratio estimate (Fig. S16e–h).

41

S7.2 DBP cost decomposition and acceleration

a

b

c

d

e

Figure S17: Additional performance analysis for the DBP application. a, Frequency-domain update time as a function of FFT size Nt for Regular and sp-DBA runs with average active-block ratios of 5%, 10%, 15%, and 20%. b, Update-stage speedup for the same Nt sweep. c, Estimated and measured end-to-end acceleration as a function of Nt . d, Normalized frequency-domain update cost as a function of the number of coupled components M for Regular and sp-DBA. e, Normalized total wall-clock cost decomposition as a function of M , separating FFT, nonlinear time-domain step, Regular frequency-domain update, and sp-DBA frequency-domain update costs.

As Nt increases, the Regular frequency-domain update grows faster than the sp-DBA update, giving larger update-stage speedups for lower active-block-ratio settings (Fig. S17a,b). Under the default 10% setting used for the main-text performance summary, the update-stage speedup reaches 21.3× and the measured end-to-end acceleration reaches 8.4× at Nt = 218 . The end-toend acceleration is bounded by unchanged FFT/IFFT calls, the nonlinear step, and sp-DBA control overhead (Fig. S17c). The M sweep shows that the savings remain concentrated in the frequency-domain linear update (Fig. S17d,e).

S8 Multi-GPU sp-DBA with distributed FFTs The distributed PFC tests apply sp-DBA to the local transformed-field update within a slabdecomposed FFT workflow. The distributed solver uses a semi-implicit pseudo-spectral update. These experiments assess distributed-workflow compatibility and update-stage scaling. Unless stated otherwise, the physical parameters, block size, activation settings, and initialization follow Supplementary Section S3.2.

42

Table S3: Core settings for the distributed PFC tests.

Item

Setting

GPU platform CUDA / MPI

NVIDIA A800-SXM4 GPUs CUDA 12.2; Open MPI 4.1.6 with CUDA-aware device-buffer communication Custom slab-decomposed 2D FFT implementation P = 1, 2, 4, 8 N = 4096 N = 2048, 2880, 4096, 5760 10,000 Three repeated runs after warm-up GPU-synchronized wall-clock time

Distributed FFT Tested GPU counts Strong-scaling grid Weak-scaling grids Time steps Repeats Timing

S8.1 Distributed FFT workflow The distributed PFC implementation combines local cuFFT real-to-complex and complex-to-real transforms with MPI slab redistribution [6, 7, 8]. Half-spectrum interpretation uses Hermitian symmetry and the remapping rule in Supplementary Section S1. Distributed workflow of sp-DBA under slab-decomposed FFT 1. Local slabs (domain partitioning)

Rank 0

...

Rank 1

...

Rank 2

...

2. Distributed forward FFT

...

4. sp-DBA active block update

5. Distributed inverse FFT

Local FFT

Local IFFT

Pack & local transpose

Pack & local transpose

MPI all-to-all

MPI all-to-all

... Rank P-1

3. Local blocks with MPI state halo

Unpack & Reassembly

Second FFT

... ...

Unpack & Reassembly

Second IFFT MPI state halo

Figure S18: Distributed workflow of sp-DBA under slab decomposition. The workflow consists of local slab partitioning, distributed forward FFT, local block construction with MPI state halos, sp-DBA transformed-field update, and distributed inverse FFT. The sp-DBA update is applied after each rank owns its assembled spectral data, while the forward and inverse distributed FFT stages retain the communication pattern required by the transform.

43

a

b

Initial domain partitioning & Row-wise 1D FFT

Pack & local transpose pack

0

1

2

3

Rank 0

B0,0

B0,1

B0,2

B0,3

x

1D FFT

Rank 0

BT0,0

BT0,1

BT0,2

BT0,3

Rank 1

B1,0

B1,1

B1,2

B1,3

x

1D FFT

Rank 1

BT1,0

BT1,1

BT1,2

BT1,3

Rank 2

B2,0

B2,1

B2,2

B2,3

x

1D FFT

Rank 2

BT2,0

BT2,1

BT2,2

BT2,3

Rank 3

B3,0

B3,1

B3,2

B3,3

x

1D FFT

Rank 3

BT3,0

BT3,1

BT3,2

BT3,3

c

d

MPI all-to-all redistribution pack

0

1

2

3

Rank 0

BT0,0

BT1,1

BT2,0

BT3,0

Rank 1

example of intra-block index remapping 1,0

0,1

1,1

Unpack & Reassembly & Column-wise 1D FFT

Rank 0

BT0,0

BT1,1

BT2,0

BT3,0

y

BT0,1

Rank 1

BT0,1

BT1,1

BT2,1

BT3,1

y

Rank 2

BT0,2

Rank 2

BT0,2

BT1,2

BT2,2

BT3,2

y

Rank 3

BT0,3

Rank 3

BT0,3

BT1,3

BT2,3

BT3,3

y

representative send paths

0,0

1D FFT 1D FFT 1D FFT 1D FFT

other similar paths

Figure S19: Slab decomposition and block redistribution in the distributed two-dimensional FFT. a, Row-wise slab decomposition before the first local FFT. b, Local block organization and transpose before MPI redistribution. c, MPI all-to-all redistribution of transposed blocks. d, Reassembled layout used for the second one-dimensional FFT.

After the distributed FFT assembles the spectral data on each rank, each rank performs the local sp-DBA update and exchanges only the block states required at rank boundaries (Fig. S19); the global transpose and all-to-all redistribution remain unchanged.

Neighbor state exchange across MPI ranks

Local boundary state update with imported halo

Active state Halo block state

State update

Dynamic block

MPI state halo

Figure S20: Rank-boundary metadata exchange for distributed sp-DBA. Block states at neighboring rank boundaries are exchanged as metadata. The imported halo states are then used to update local boundary states so that active and transition blocks are assigned consistently across rank partitions.

The distributed implementation uses tiled slab packing, self-rank bypass, CUDA-aware 44

nonblocking MPI, and device-buffer metadata exchange, while the sp-DBA update mechanism itself is unchanged.

S8.2 Strong and weak scaling

a

b

Strong scaling

c

Weak scaling

d

Figure S21: Distributed scaling of sp-DBA in the PFC application. a, Update stage under strong scaling at fixed global grid size N = 4096. b, Update stage under weak scaling with approximately fixed local grid area per GPU. c, Overall strong-scaling wall-clock decomposition at fixed global grid size N = 4096. d, Overall weak-scaling wall-clock decomposition with approximately fixed local grid area per GPU. Gray and blue bars denote the Regular and sp-DBA update stages, respectively; numbers above the blue bars indicate the measured average active-block ratio, not speedup.

Under strong scaling, both update times decrease with the local workload, with the sp-DBA update time below the Regular update time (Fig. S21a). Under weak scaling, the sp-DBA update time changes little and the measured average active-block ratio stays close to 10% (Fig. S21b). Distributed FFT and communication dominate the total wall-clock cost (Fig. S21c,d), so the overall acceleration is bounded by the share of the local update in the distributed workflow.

45

S8.3 CUDA-aware MPI communication

a

b

Figure S22: Communication path comparison between host staging and CUDA-aware MPI. a, Strongscaling control case at P = 8 and N = 4096. b, Weak-scaling control case at P = 8 and N = 5760. Within each solver mode, the host-staged path is shown on the left and CUDA-aware device-buffer MPI is shown on the right. Numbers above the bars indicate total wall-clock time.

Reported timings use CUDA-aware device-buffer MPI; host-staged runs serve only as communicationpath controls (Fig. S22).

S8.4 Activation stability under distributed execution

a

b

c

Figure S23: Activation stability in distributed PFC tests. a, Measured active-block ratio as a function of the target average active-block ratio at P = 1 and P = 8. The dashed line denotes exact target tracking. b, Large-N update-stage comparison at N = 8192 on 4 and 8 GPUs. Gray and blue bars denote the Regular and sp-DBA update times, respectively; numbers above the blue bars indicate the measured average active-block ratio. c, Average active-block ratio during a 50,000-step distributed sp-DBA run at P = 8 and N = 4096.

The target-ratio test varies the target average active-block ratio from 5% to 20%. The measured ratio follows the target but remains below the exact tracking line in the distributed setting (Fig. S23a). The 10% setting is used as the default because it provides a clear update-time reduction while keeping the measured ratio near the target range. At N = 8192, sp-DBA reduces the local transformed-field update time on both 4 and 8 GPUs relative to Regular (Fig. S23b). Over 50,000 steps at P = 8 and N = 4096, the active-block ratio stays close to the default 46

10% range, indicating stable dynamic threshold control and rank-boundary state exchange under the tested decomposition (Fig. S23c).

S9 Visualization ParaView [75] was used for field, Fourier-amplitude, isosurface, and volume rendering, while OVITO [76] was used for XPFC grain reconstruction and morphology visualization. Activation maps were generated from recorded block labels, while scalar quantities and performance plots were prepared from exported analysis files. Post-processing was excluded from solver timing unless stated otherwise.

47

Related documents

• sp-DBA: a general framework for adaptive transform-domain computation#926788
• Understanding Domain-Aware Distribution Alignment in Budgeted Entity Matching#310799
• From IOCs to Regex: Automating CTI Operationalization for SOC with LLMs#13029
• DataFlow: An LLM-Driven Framework for Unified Data Preparation and Workflow Automation in the Era of Data-Centric AI#192040
• jaxDecomp: JAX Library for 3D Domain Decomposition and Parallel FFTs#325939