Communication Strategy Selection for Multi-GPU 3D FDTD with Convolutional Perfectly Matched Boundary Layers
arXiv:2606.06910v1 [cs.DC] 5 Jun 2026
Victory Obieke Department of Mathematics Oregon State University Corvallis, OR, USA [email protected]
Abstract
per GPU is small [2, 5, 6]. Enlarged ghost regions reduce the communication frequency by exchanging a wider halo and advancing several local time steps before the next exchange. This idea is related to temporal blocking and ghost-zone optimization for stencil computations [2, 5, 13]. Most idealized stencil benchmarks omit boundary treatments that are essential in production FDTD solvers. Perfectly matched layers and convolutional perfectly matched layers (CPML) are widely used to reduce artificial reflections at computational boundaries [1, 10, 11]. CPML changes the performance balance because it introduces auxiliary variables, boundary-layer updates, and additional memory traffic. Therefore, the benefit of enlarged ghost regions must be evaluated in the presence of realistic absorbing boundary costs. The goal of this work is not to show that four GPUs are always faster than one GPU. Instead, the goal is to determine which communication strategy is most effective when a 3D FDTD+CPML solver is distributed across multiple GPUs. The contribution of this work is an empirical communicationstrategy study showing that, for high-order 3D FDTD+CPML on peer-connected GPUs, direct GPU-to-GPU exchange is the dominant optimization, while enlarged ghost regions provide only limited additional benefit. We first compare baseline decomposition layouts, and then use the selected pencil-yz layout to study host-staged exchange, direct GPU-to-GPU peer exchange, and enlarged ghost-region communication. The results show that direct GPU-to-GPU peer exchange is the most important optimization, enlarged ghost regions give only modest additional speedup, and multi-GPU decomposition is most valuable when larger grid sizes approach or exceed the memory capacity of a single GPU.
In this paper we describe a communication-strategy study for multi-GPU three-dimensional finite-difference time-domain computation with convolutional perfectly matched layer boundary conditions using CUDA. The metrics used to determine the most effective implementation include runtime, throughput in millions of output points per second, strongscaling efficiency, CPML overhead, host-staged versus direct GPU-to-GPU exchange speedup, and enlarged-ghost speedup. On a single NVIDIA Quadro RTX 6000 GPU, the CPML implementation sustains 2,889–3,290 million output points per second with less than 1% boundary-layer overhead, providing the single-GPU baseline for the multi-GPU study. The results show that direct GPU-to-GPU peer exchange is the dominant optimization with a 2.46–2.76× speedup over hoststaged exchange, while enlarged ghost regions give only modest benefits because the reduced communication frequency is partly offset by redundant computation and additional memory traffic. On NVIDIA Quadro RTX 8000 GPUs, the implementation gives up to a 1.51× speedup on two GPUs for the tested strong-scaling cases, while four GPUs enable larger grids that approach or exceed single-GPU memory capacity.
1
Introduction
Finite-difference time-domain (FDTD) methods are widely used for wave propagation, electromagnetics, seismic modeling, and computational physics because they combine structured-grid simplicity with regular stencil updates [4, 11, 12]. This structure is well suited for GPUs, where many grid points can be updated in parallel and large memory bandwidth can be exploited [3,6,7]. However, practical three-dimensional simulations often require large domains, high-order stencils, and absorbing boundary layers, so the memory footprint can exceed the capacity of a single GPU [6, 10, 11]. A central challenge in multi-GPU stencil computation is the balance between local computation and inter-device communication. In the standard one-step halo exchange, neighboring GPUs exchange ghost layers after every time step. This method is simple and avoids redundant work, but it can become communication dominated when the local subdomain
2
Related Work
GPU acceleration of high-order three-dimensional finitedifference stencils has been studied extensively. Micikevicius [6] used data-access redundancy as a key metric for optimizing 3D finite-difference computation on CUDA GPUs and also discussed multi-GPU extensions. Communicationavoiding and temporal-blocking methods reduce communication or memory traffic by computing multiple time steps lo1
3.2
cally, but they introduce redundant work in overlap or ghost regions [2, 5, 13]. The CPML formulation of Roden and Gedney provides an efficient implementation of complex-frequencyshifted PML for FDTD simulations [10]. In contrast to interior-only stencil benchmarks, the present work includes CPML boundary layers and measures their effect in multiGPU runs. Table 1 summarizes the positioning of this work relative to closely related FDTD, CPML, and GPU stencilcomputation studies. The novelty of the present work is not the introduction of a new FDTD scheme, but the combined performance evaluation of decomposition layout, host-staged communication, direct GPU-to-GPU peer exchange, and enlarged ghost-region communication for a practical multi-GPU 3D FDTD solver with CPML boundary layers.
3
Model
3.1
First-Order Acoustic Model and Discrete Update
CPML Boundary Layers
The absorbing boundary layer follows the CFS/Roden– Gedney-style CPML construction [10]. In the interior of the domain, the first-order acoustic FDTD system is updated using the usual finite-difference derivatives. Inside the CPML layer, however, each spatial derivative is replaced by a stretched derivative plus a recursive memory correction. Thus, a discrete derivative such as δx f is replaced by (1/κx )δx f + ψx , where f can be the pressure p or one of the velocity components vx , vy , and vz . The variable ψx is the CPML memory variable associated with the x-direction derivative. The memory variable is updated recursively by ψxn+1 = bx ψxn + ax δx f n . This recursive update is the practical implementation of the convolutional correction in CPML. In other words, the CPML correction depends on the history of the spatial derivative, but the full history is not stored. Instead, the single memory variable ψx carries the required history information efficiently from one time step to the next. The CPML coefficients are computed from damping and stretching profiles inside the absorbing layer. Let ρ ∈ [0, 1] denote the normalized distance into the CPML layer. In the performance code, cubic grading is used, so q = 3 and ρq = ρ3 . The damping profile is σx (ρ) = σmax ρ3 , where σmax = −4 log(R)/(2LCPML h) and R = 10−8 . The stretching profile is κx (ρ) = 1 + (κmax − 1)ρ3 , with κmax = 5. The complex-frequency-shift profile is αx (ρ) = αmax (1 − ρ), where αmax = 0.05σmax . Analogous profiles are used in the y- and z-directions. The recursive coefficients are bx = exp[−(σx /κx + αx )∆t] and ax = σx (bx −1)/[κx (σx +κx αx )]. The coefficient bx controls the decay of the memory variable, while ax controls how strongly the current derivative δx f n contributes to the memory variable. Outside the CPML layer, σx = 0, κx = 1, and the memory variable is zero, so the modified derivative (1/κx )δx f + ψx reduces to the ordinary finite-difference derivative δx f . Therefore, the CPML modification does not change the interior FDTD update. For the velocity update, the standard interior update n−1/2 n+1/2 − ∆t δx pn is replaced inside the CPML vx = vx n+1/2 n−1/2 n by vx = vx − ∆t[(1/κx )δx pn + ψv,x ]. Similarly,
In the interior of the domain, away from the absorbing layers, we use the first-order acoustic pressure–velocity system pt = −K∇ · v and vt = −(1/ρ)∇p, where p is pressure, v = (vx , vy , vz ) is particle velocity, K is the bulk modulus, and ρ is density. This system implies the scalar wave equation ptt = c2 ∆p, where c2 = K/ρ. For the performance experiments in this work, the coefficients are absorbed into the wave-update scaling parameter, so the update is written in normalized form. (2r) (2r) (2r) Let δx , δy , and δz denote centered finite-difference derivative operators of spatial order 2r. The spatial derivatives in the acoustic system are approximated by centered (2r) (2r) (2r) finite-difference operators. We write δx , δy , and δz for the numerical approximations of ∂x , ∂y , and ∂z , respec(2r) tively. The Pr x-derivative is approximated by (δx f )i,j,k = (1/∆x) ℓ=1 aℓ (fi+ℓ,j,k − fi−ℓ,j,k ). Similarly, the y(2r) and z-derivatives are approximated by (δy f )i,j,k = Pr (2r) (1/∆y) Pℓ=1 aℓ (fi,j+ℓ,k − fi,j−ℓ,k ) and (δz f )i,j,k = r (1/∆z) ℓ=1 aℓ (fi,j,k+ℓ − fi,j,k−ℓ ). Here, r is the stencil radius and aℓ are the finite-difference weights. In the eighthorder experiments, r = 4. Without CPML terms, the staggered pressure–velocity upn+1/2 n−1/2 date is written componentwise as vx = vx − (2r) n n+1/2 n−1/2 (2r) n (∆t/ρ)δx p , vy = vy − (∆t/ρ)δy p , and n+1/2 n−1/2 (2r) vz = vz − (∆t/ρ)δz pn . The pressure is then up(2r) n+1/2 (2r) n+1/2 dated by pn+1 = pn − K∆t(δx vx + δ y vy + (2r) n+1/2 δ z vz ). Thus one full time step consists of a pressure-to-velocity stencil update followed by a velocity-to-pressure stencil update. Since each subupdate uses a radius-r stencil, numerical dependence can extend by 2r grid cells during one full time step. Therefore, if s local time steps are taken between halo exchanges, the enlarged ghost depth is chosen as g = 2rs. This choice provides enough neighboring data for the local subdomain to advance s steps before the next communication event.
n+1/2 n−1/2 n+1/2 n = vy − ∆t[(1/κy )δy pn + ψv,y ], and vz = n−1/2 n n vz − ∆t[(1/κz )δz p + ψv,z ]. Here, ψv,x , ψv,y , and ψv,z
vy
are the CPML memory variables associated with the pressure derivatives used in the velocity updates. The pressure update is modified in the same way. In the interior, pressure is updated from the divergence of the velocity field. Inside the CPML layer, this update becomes n+1/2 n+1/2 pn+1 = pn − ∆t[(1/κx )δx vx + (1/κy )δy vy + n+1/2 n+1/2 n+1/2 n+1/2 (1/κz )δz vz + ψp,x + ψp,y + ψp,z ]. The memory variables ψp,x , ψp,y , and ψp,z correspond to the velocity derivatives used in the pressure update. Thus, each coordinate direction has its own damping profile, stretching coefficient, and memory correction. The CPML update is applied only in the global absorbing boundary layer. At internal GPU interfaces, no CPML 2
Table 1: Positioning of the present work relative to closely related GPU stencil and FDTD studies. Study or research direction
Multi-GPU CPML/PML High-order 3D FDTD Decomposition comparison Host vs. peer exchange Enlarged ghost regions
Classical FDTD and PML formulations [1, 10, 12]
✗
✓
✗
✗
✗
✗
GPU stencil optimization and temporal blocking studies [2, 5, 13]
✓
✗
✓
✗
✗
✓
GPU 3D finite-difference / FDTD acceleration studies [6, 7]
✓
✗
✓
✗
✗
✗
CPML-based practical FDTD simulations [10, 11]
✗
✓
✓
✗
✗
✗
Present work
✓
✓
✓
✓
✓
✓
is used. Neighboring GPU subdomains communicate only through halo exchange, while the CPML memory variables are updated only near the physical outer boundary of the global computational domain. This separation is important because CPML represents an absorbing physical boundary treatment, whereas halo exchange is only a communication mechanism used to provide neighboring grid values across internal GPU subdomain interfaces.
3.3
Multi-GPU Decomposition and Halo Exchange
The three-dimensional grid is decomposed across four GPUs using three layouts: slab-z, block-xy, and pencil-yz, corresponding to 1 × 1 × 4, 2 × 2 × 1, and 1 × 2 × 2, respectively. In the slab decomposition, the domain is split only in the zdirection. In the block decomposition, the domain is split in x and y. In the pencil-yz decomposition, the domain is split in y and z, while the full x-direction remains local to each GPU. Each GPU stores its owned subdomain plus ghost cells in the decomposed directions. Physical outer boundaries are treated by CPML, while internal GPU interfaces are treated only by halo exchange. Thus, CPML is applied only at the global outer boundary of the computational domain, not at internal subdomain interfaces. The implementation uses raw CUDA kernels through CuPy to pack halo data into contiguous send buffers, direct CUDA peer copies for GPU-to-GPU exchange when available, and raw CUDA kernels to unpack the received data into ghost cells. The local FDTD+CPML update is then applied on each GPU. For standard halo exchange, s = 1, ghost regions are exchanged every time step. For enlarged ghost exchange, each GPU exchanges a deeper halo and then advances s local steps before the next exchange. Since one full first-order acoustic time step contains both a pressure-to-velocity update and a velocity-to-pressure update, the enlarged ghost depth is chosen as g = 2rs. This reduces communication frequency but increases redundant computation and memory traffic in the enlarged ghost region.
Figure 1: Four-GPU pencil-yz decomposition.
performance experiments use single precision to emphasize communication and memory-bandwidth behavior. Reported runtimes are medians over repeated runs.
4.1
Hardware and Software Environment
The reported GPU experiments were run on the optimus node in the preempt Slurm partition. The node contains two Intel Xeon Gold 6230 CPUs and four NVIDIA Quadro RTX 6000 GPUs. Because communication results depend on the node topology, we record the CUDA software stack and the GPU-to-GPU peer-access topology. In this topology, PIX denotes communication through at most one PCIe bridge, NODE denotes communication through PCIe and host bridges within a NUMA node, and SYS denotes communication that traverses PCIe and the inter-socket CPU interconnect. The peer-exchange results should therefore be interpreted as CUDA peer-access results on a PCIe/NUMAconnected four-GPU node, not as NVLink results. Before applying the final raw-kernel communication optimizations, we first compare the three four-GPU decomposition layouts. This baseline comparison is used only to select the most favorable domain decomposition for the subsequent 4 Experimental Setup optimized communication experiments. Experiments use a raw-kernel implementation of the firstTable 5 shows that the pencil-yz decomposition gives the order acoustic FDTD+CPML system. Unless otherwise stated, best throughput for all tested grids in the baseline implemen3
Table 2: Numerical, CPML, and communication parameters used in the multi-GPU performance experiments.
Table 4: CUDA peer-access matrix and topology labels for the selected four-GPU node. A value of 1 indicates that CUDA peer access is available. Topology labels are taken from nvidia-smi topo -m.
Parameter Value Governing system First-order acoustic pressure–velocity FDTD system Spatial stencil order Eighth order Stencil radius r=4 Time discretization Two-stage pressure–velocity update Wave-update scaling / CFL parameter λFDTD = 0.05 Communication intervals s ∈ {1, 2, 4, 8} Enlarged ghost depth g = 2rs Boundary condition CFS/Roden–Gedney-style CPML CPML application region Global outer boundary only Internal GPU interfaces Halo exchange only, no CPML CPML thickness LCPML = 20 cells CPML target reflection parameter R = 10−8 CPML polynomial grading Cubic grading, q = 3 CPML damping profile σ(ρ) = σmax ρ3 CPML maximum damping σmax = −4 log(R)/(2LCPML h) CPML stretching profile κ(ρ) = 1 + (κmax − 1)ρ3 CPML maximum stretching parameter κmax = 5 CPML complex-frequency shift α(ρ) = αmax (1 − ρ) CPML shift parameter αmax = 0.05σmax CPML auxiliary fields Recursive memory variables for pressure and velocity derivatives Performance precision Single precision Verification precision Double precision where stated Runtime statistic Median over repeated runs
From/To GPU 0 GPU 1 GPU 2 GPU 0 – 1 / NODE 1 / NODE GPU 1 1 / NODE – 1 / PIX GPU 2 1 / NODE 1 / PIX – GPU 3 1 / SYS 1 / SYS 1 / SYS
Table 5: Baseline four-GPU decomposition comparison for the first-order acoustic 3D FDTD+CPML benchmark. This table uses the pre-optimized communication implementation. Grid slab-z (1 × 1 × 4) block-xy (2 × 2 × 1) pencil-yz (1 × 2 × 2) 3203 2.357 s / 889.71 3.417 s / 613.72 2.345 s / 894.45 5443 8.399 s / 1226.79 10.889 s / 946.24 7.616 s / 1352.80 8003 23.778 s / 1378.06 29.947 s / 1094.19 21.932 s / 1494.09
Table 3: Hardware and software environment for the main fourGPU RTX 6000 experiments.
decompositions. The reduction happens with more layers as expected in [10].
Component Value Node optimus.hpc.engr.oregonstate.edu Slurm partition preempt CPU 2× Intel Xeon Gold 6230, 20 cores/socket CPU cores 40 physical cores, 2 NUMA nodes System memory 754 GiB GPU model NVIDIA Quadro RTX 6000 Number of GPUs 4 GPU memory 23040 MiB nominal; 22501 MiB available to CUDA per GPU GPU driver 590.48.01 CUDA version reported by nvidia-smi 13.1 CUDA toolkit 13.0, V13.0.88 CUDA runtime / driver API 13020 / 13010 CuPy version 14.1.0 Python version 3.10.14 Operating system Linux 5.14.0-570.58.1.el9 6.x86 64 GPU topology PCIe/NUMA paths with PIX, NODE, and SYS links; no NVLink reported Peer access Enabled for all ordered GPU pairs Communication backend CuPy CUDA copies with peer access enabled NCCL Not used
To complement the reflection-ratio verification, we also include a qualitative wavefield visualization of the absorbing behavior of CPML. This experiment follows the acoustic CPML validation style of Pasalic and McGarry [9], who derive CPML for isotropic and anisotropic acoustic wave equations and demonstrate its effectiveness using source-driven acoustic wave propagation. The purpose of the present visualization is not to reproduce their full TTI anisotropic model exactly, but to provide a comparable acoustic CPML demonstration for the pressure–velocity FDTD setting used in this work. The visualization uses a two-dimensional acoustic-style pressure–velocity system on a 2 km × 2 km domain with 2 ∆x = ∆z = 8 m. Source s(t) = (1 − 2a2 )e−a , with a = πf0 (t − t0 ), is placed at the center of the domain with peak frequency f0 = 30 Hz. We set cx = 2600 m/s and cz = 1800 m/s. The CPML thickness is 10 grid cells. The same source, grid, time step, and final time are used for the no-CPML and CPML simulations.
tation. Therefore, the pencil-yz layout shown in Figure 1 is used for the remaining optimized four-GPU experiments.
5
GPU 3 1 / SYS 1 / SYS 1 / SYS –
FDTD and CPML Implementation Verification
Figure 2 compares pressure snapshots with and without CPML. In the no-CPML case, waves reflect from the artificial To verify the finite-difference implementation, we boundary and remain visible inside the computational domain. use the periodic standing wave u(x, = In the CPML case, the outgoing wave is strongly attenuated as √ y, z, t) sin(2πmx) sin(2πmy) sin(2πmz) cos(2πm 3 ct), which it enters the absorbing layer, leaving a much smaller reflected satisfies utt = c2 ∆u using the method of manufactured field. Figure 3 shows the corresponding normalized field energy. After the source injection ends, the no-CPML simulation solutions Tables 6 and 7 show the expected spatial and temporal or- retains energy because waves are trapped by boundary reflections, while the CPML simulation decays by several orders of ders of accuracy. magnitude as outgoing waves are absorbed.
5.1
FDTD Verification
Figures 2 and 3 show a separate single-domain acoustic CPML demonstration. They are included only to illustrate the We further verify the CPML implementation in Table 8. absorbing behavior of CPML visually, while Table 8 provides Table 8 confirms that the CPML reduces artificial reflections the decomposition-dependent quantitative check used for the by several orders of magnitude and behaves consistently across multi-GPU implementation.
5.2
CPML Verification
4
Table 6: Spatial convergence of the four-GPU pencil-yz finitedifference Laplacian for the periodic standing-wave solution with mode number m = 4. Stencil order
2
4
8
N Relative L2 error 30 5.714 × 10−2 60 1.454 × 10−2 120 3.650 × 10−3 240 9.135 × 10−4 480 2.284 × 10−4 30 5.141 × 10−3 60 3.368 × 10−4 120 2.130 × 10−5 240 1.335 × 10−6 480 8.349 × 10−8 30 6.568 × 10−5 60 2.891 × 10−7 120 1.164 × 10−9 240 4.561 × 10−12
Rate Relative L∞ error – 5.714 × 10−2 1.97 1.454 × 10−2 1.99 3.650 × 10−3 2.00 9.135 × 10−4 2.00 2.284 × 10−4 – 5.141 × 10−3 3.93 3.368 × 10−4 3.98 2.130 × 10−5 4.00 1.335 × 10−6 4.00 8.349 × 10−8 – 6.568 × 10−5 7.83 2.891 × 10−7 7.96 1.164 × 10−9 7.99 4.606 × 10−12
Table 8: Reflected-to-incident amplitude ratio for different fourGPU decompositions. CPML cells 0 8 12 16 24
Rate – 1.97 1.99 2.00 2.00 – 3.93 3.98 4.00 4.00 – 7.83 7.96 7.98
slab-z (1 × 1 × 4) 9.438 × 10−1 8.049 × 10−2 6.714 × 10−3 5.581 × 10−4 3.328 × 10−5
block-xy (2 × 2 × 1) 9.438 × 10−1 8.049 × 10−2 6.715 × 10−3 5.581 × 10−4 5.300 × 10−5
pencil-yz (1 × 2 × 2) 9.438 × 10−1 8.049 × 10−2 6.714 × 10−3 5.581 × 10−4 3.328 × 10−5
Table 7: Time convergence of the four-GPU pencil-yz wave solver using the eighth-order spatial stencil and mode number (m=1). ∆t 1.0 × 10−3 5.0 × 10−4 2.5 × 10−4
6
Relative L2 error 1.559 × 10−5 3.887 × 10−6 9.704 × 10−7
Rate – 2.00 2.00
Relative L∞ error 1.559 × 10−5 3.887 × 10−6 9.704 × 10−7
Rate – 2.00 2.00
Performance Metrics
Runtime is the wall-clock time required to advance the full 3D grid for Nt time steps. Throughput is reported as Mpoints/s = Nx Ny Nz Nt /(106 Truntime ). Strong-scaling efficiency on p GPUs is Ep = T1 /(pTp ), and the enlarged-ghost Figure 2: Acoustic CPML wavefield visualization using a Ricker speedup relative to s = 1 is Ss = Ts=1 /Ts . CPML overhead pulse source. The no-CPML case retains reflected waves, while the CPML case strongly attenuates the outgoing wave and reis computed as (TCPML − Tno CPML )/Tno CPML × 100%. duces boundary reflections.
7
Results and Discussion
7.1
Decomposition Comparison
For the multi-GPU runs, each GPU stores an owned subdomain together with ghost cells along internal GPU interfaces. Before each local update, halo data are packed into contiguous device buffers using raw CUDA packing kernels. These buffers are then transferred directly between neighboring GPUs using CUDA peer-to-peer copies, and unpacked into the receiving ghost cells using raw CUDA unpacking kernels. This avoids host-staged communication of the form GPU–CPU–GPU and instead uses direct GPU–GPU exchange whenever peer access is available. The optimized code also reuses communication buffers rather than allocating new buffers at every time step. For a standard halo exchange, s = 1, the halo is exchanged every time step. For enlarged ghost-region communication, each GPU exchanges a deeper ghost region of depth g = 2rs, where r is the stencil radius and s is the number of local steps taken before the next exchange. This reduces the number of halo exchanges by a factor of approximately s, but it also introduces redundant computation in the enlarged ghost region. The performance tables below use this raw-kernel communication path for the pencil-yz layout. With the decomposition fixed, the experiments isolate the effect of communica-
Most experiments were performed on the four-GPU NVIDIA Quadro RTX 6000 node described in Table 3. The strongscaling experiment in Table 11 was performed separately on NVIDIA Quadro RTX 8000 GPUs, and is reported separately because it uses different hardware. Each entry gives runtime in seconds followed by throughput in Mpoints/s.
7.2
Raw-Kernel and Communication Optimizations
The optimized implementation used in the performance experiments replaces high-level array operations with explicit raw CUDA kernels launched through CuPy. Each kernel assigns CUDA threads directly to grid points in the local subdomain, so the pressure field, velocity fields, and CPML auxiliary memory variables are updated using contiguous device arrays. This reduces Python-level overhead and avoids repeated temporary-array creation during the time-stepping loop. 5
Table 9: Correctness check for enlarged ghost-region communication in the four-GPU pencil-yz first-order acoustic FDTD+CPML implementation. For each grid, the table reports the maximum relative difference over s ∈ {2, 4, 8} compared with the standard s = 1 halo-exchange reference solution. Grid 3203 4803 5443 6403 8003
GPUs 4 4 4 4 4
Max relative L2 difference 0.000 × 100 0.000 × 100 0.000 × 100 0.000 × 100 0.000 × 100
Max relative L∞ difference 0.000 × 100 0.000 × 100 0.000 × 100 0.000 × 100 0.000 × 100
Table 10: Single-GPU raw-kernel throughput for the firstorder acoustic 3D FDTD benchmark with and without CPML boundary-layer updates.
Figure 3: Normalized field energy with and without CPML. After the source injection ends, the no-CPML case retains energy due to boundary reflections, while the CPML case decays by several orders of magnitude as outgoing waves are absorbed.
Grid Boundary Runtime (s) Mpoints/s CPML overhead 3203 none 0.635 3304.48 – 3203 CPML 0.638 3289.58 0.45% 4803 none 2.440 2900.31 – 4803 CPML 2.449 2889.58 0.37% 5443 none 3.530 2918.65 – 5443 CPML 3.553 2900.14 0.64% 6403 none 5.751 2917.52 – 3 640 CPML 5.766 2909.64 0.27% 8003 none 11.204 2924.73 – 8003 CPML 11.221 2920.36 0.15%
are expected because the overhead is computed from small differences between separate runtime measurements.
7.5
Table 11 shows two regimes. For grids up to 8003 , the best runtime is obtained with two GPUs; four GPUs still gives speedup over one GPU for the 5443 and 8003 cases, but it is slower than the two-GPU configuration because the additional halo-exchange and synchronization costs outweigh the extra reduction in local work. For larger grids, the memorycapacity benefit becomes more important. The 10243 case runs slightly faster on four GPUs than on two GPUs. Thus, additional GPUs become useful not only for runtime reduction, but also for enabling larger CPML simulations that do not fit on fewer devices.
Figure 4: Halo exchange workflow.
tion strategy: standard one-step halo exchange, enlarged ghost regions, direct GPU-to-GPU peer exchange, and host-staged exchange.
7.3
Correctness Check
Table 9 shows that the enlarged ghost-region solutions agree with the standard s = 1 reference to the reported precision. The zero values in Table 9 mean that the differences are below the printed precision; in these tests the enlarged-ghost solutions were numerically indistinguishable from the standard s = 1 reference.
7.4
Strong Scaling and Memory Capacity on RTX 8000 GPUs
7.6
Enlarged Ghost-Region Performance
After selecting the pencil-yz layout from the baseline decomposition study in Table 5, we apply the optimized raw-kernel communication implementation and study the effects of enlarged ghost regions and direct GPU-to-GPU peer exchange. Table 12 shows that enlarged ghost exchange improves performance for all tested grids,with the best performance at s = 4. The s = 8 case is slower than s = 4, indicating that excessive ghost enlargement introduces enough redundant work to reduce the benefit of fewer exchanges.
Single-GPU Baseline
Table 10 shows that the measured single-GPU CPML boundary-layer overhead remains below 1% for all tested grids. Since the CPML update is confined to a thin boundary layer, its relative cost is small compared with the full-volume FDTD update. Small nonmonotone variations across grid sizes 6
Table 11: Strong-scaling and memory-capacity behavior of the raw-kernel first-order acoustic FDTD+CPML benchmark on NVIDIA Quadro RTX 8000 GPUs using standard s = 1 halo exchange.
Table 12: Four-GPU enlarged ghost-region performance for the optimized raw-kernel pencil-yz first-order acoustic FDTD+CPML benchmark. The case s = 1 is the standard onestep halo-exchange method.
Grid
Grid GPUs 4 4 3203 4 4 4 4 3 480 4 4 4 4 5443 4 4 4 4 3 640 4 4 4 4 8003 4 4
3203
5443
8003
10243
7.7
GPUs Runtime (s) Speedup Parallel efficiency 1 0.897 1.00 100.0% 2 0.864 1.04 51.9% 4 1.600 0.56 14.0% 1 4.666 1.00 100.0% 2 3.442 1.36 67.8% 4 4.482 1.04 26.0% 1 14.614 1.00 100.0% 2 9.658 1.51 75.7% 4 10.458 1.40 34.9% 1 OOM – – 2 19.050 – – 4 18.528 – –
Host-Staged Versus Peer Exchange
Table 13 shows that direct GPU-to-GPU peer exchange gives 2.46×–2.76× speedup over host-staged exchange.
7.8
Memory Capacity Benefit
s Runtime (s) Mpoints/s Speedup vs. s = 1 1 1.708 1227.87 1.00 2 1.556 1348.07 1.10 4 1.491 1406.48 1.15 8 1.568 1337.26 1.09 1 3.669 1929.27 1.00 2 3.521 2010.15 1.04 4 3.448 2052.46 1.06 8 3.636 1946.61 1.01 1 4.791 2150.36 1.00 2 4.599 2240.38 1.04 4 4.514 2282.74 1.06 8 4.727 2179.67 1.01 1 6.781 2474.16 1.00 2 6.537 2566.60 1.04 4 6.397 2622.52 1.06 8 6.654 2521.52 1.02 1 11.100 2951.97 1.00 2 10.738 3051.65 1.03 4 10.444 3137.44 1.06 8 10.799 3034.24 1.03
Table 14 reports an analytical estimate of the per-GPU field storage. For grids that fit comfortably on one GPU, single- such as energy-conserving data-driven convolution stencils for GPU execution can be more cost-effective because it avoids Maxwell-type systems [8]. communication overhead. For larger CPML simulations, however, multi-GPU decomposition reduces the per-GPU memory Availability of Data and Code footprint and enables runs that are close to or beyond the capacity of one device. The code and scripts used to generate the numerical and performance results in this study are available at https://github.com/victoryobieke/ 8 Conclusions and Future Work fdtd-cpml-multigpu. Overall, the experiments show that the main benefit of the multi-GPU implementation is not universal strong-scaling Conflict of Interest speedup over a highly optimized single-GPU solver. Rather, the main benefits are communication efficiency and memory The author declares that there are no conflicts of interest. scalability. Direct GPU-to-GPU peer exchange removes most of the host-staging bottleneck, while multi-GPU decomposition reduces the per-GPU memory footprint and enables larger References CPML simulations. Enlarged ghost regions can further reduce communication frequency, but their benefit is limited by re- [1] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. Journal of computational dundant computation and extra memory traffic. physics, 114(2):185–200, 1994. Future work will extend the implementation to multi-node GPU systems with improved communication overlap using [2] S. Kamil, K. Datta, S. Williams, L. Oliker, J. Shalf, and CUDA streams, NCCL, or GPU-aware MPI. A second direcK. Yelick. Implicit and explicit optimizations for stencil tion is to extend the study to full Maxwell systems, heterogecomputations. In Proceedings of the 2006 workshop on neous media, higher-order CPML models, and other boundMemory system performance and correctness, pages 51– ary conditions, with the goal of automatically selecting de60, 2006. composition layouts and ghost-region depths for a given hardware topology and problem size. Another direction is to [3] E. Lindholm, J. Nickolls, S. Oberman, and J. Montrym. combine the present multi-GPU FDTD+CPML implementaNvidia tesla: A unified graphics and computing architection with learned structure-preserving spatial discretizations, ture. IEEE micro, 28(2):39–55, 2008. 7
Table 13: Comparison of host-staged and direct GPU-to-GPU peer exchange for the four-GPU raw-kernel pencil-yz first-order acoustic FDTD+CPML benchmark. All runs use the standard s = 1 halo-exchange method. Grid GPUs 4 3203 4 4 4803 4 4 3 544 4 4 3 640 4 4 8003 4
trary media. Microwave and optical technology letters, 27(5):334–339, 2000. [11] A. Taflove, S. C. Hagness, and M. Piket-May. Computational electromagnetics: the finite-difference timedomain method. The Electrical Engineering Handbook, 3(629-670):15, 2005.
Exchange method Runtime (s) Mpoints/s Speedup Host-staged 4.720 444.33 1.00 GPU-to-GPU peer 1.708 1227.87 2.76 Host-staged 9.628 735.17 1.00 GPU-to-GPU peer 3.669 1929.27 2.62 Host-staged 12.223 842.97 1.00 GPU-to-GPU peer 4.791 2150.36 2.55 Host-staged 16.853 995.48 1.00 GPU-to-GPU peer 6.781 2474.16 2.49 Host-staged 27.278 1201.24 1.00 GPU-to-GPU peer 11.100 2951.97 2.46
[12] K. Yee. Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media. IEEE Transactions on antennas and propagation, 14(3):302–307, 1966. [13] L. Zhang, M. Wahib, P. Chen, J. Meng, X. Wang, T. Endo, and S. Matsuoka. Revisiting temporal blocking stencil optimizations. In Proceedings of the 37th International Conference on Supercomputing, pages 251–263, 2023.
Table 14: Estimated per-GPU memory footprint for the rawkernel first-order acoustic FDTD+CPML implementation using standard s = 1 halo exchange. The estimate includes the pressure field, three velocity fields, and six CPML auxiliary memory fields in single precision. Grid GPUs Decomposition Estimated memory/GPU Memory reduction 1 1×1×1 1.41 GiB 1.00× 3203 2 1×1×2 0.74 GiB 1.91× 4 1×2×2 0.39 GiB 3.64× 1 1×1×1 6.54 GiB 1.00× 5443 2 1×1×2 3.36 GiB 1.94× 4 1×2×2 1.73 GiB 3.78× 1 1×1×1 20.24 GiB 1.00× 8003 2 1×1×2 10.32 GiB 1.96× 4 1×2×2 5.26 GiB 3.85×
[4] G. A. McMechan. Migration by extrapolation of timedependent boundary values. 2008. [5] J. Meng and K. Skadron. Performance modeling and automatic ghost zone optimization for iterative stencil loops on gpus. In Proceedings of the 23rd international conference on Supercomputing, pages 256–265, 2009. [6] P. Micikevicius. 3d finite difference computation on gpus using cuda. In Proceedings of 2nd workshop on general purpose processing on graphics processing units, pages 79–84, 2009. [7] J. Nickolls, I. Buck, M. Garland, and K. Skadron. Scalable parallel programming with cuda: Is cuda the parallel programming model that application developers have been waiting for? Queue, 6(2):40–53, 2008. [8] V. Obieke and A. E. Sunday. An energy stable approach for learning derivative operators from noisy data for maxwells equations, 2026. [9] D. Pasalic and R. McGarry. Convolutional perfectly matched layer for isotropic and anisotropic acoustic wave equations. In SEG International Exposition and Annual Meeting, pages SEG–2010. SEG, 2010. [10] J. A. Roden and S. D. Gedney. Convolution pml (cpml): An efficient fdtd implementation of the cfs–pml for arbi8