arXiv:2609.24294v1 [cs.DC] 21 Sep 2026
Toward GPU-Resident Climate Models: A Feasibility Study on Lossy Compression for the Spherical Harmonic Transform’s Communication Bottleneck Lorenzo Breschi
Flavio Vella
University of Trento Trento, Italy [email protected]
University of Trento Trento, Italy [email protected]
Abstract—Operational pseudospectral atmospheric models such as the ECMWF Integrated Forecasting System (IFS) run today almost exclusively on CPUs; GPU ports are under active development but not yet used in production. These models rely on the Spherical Harmonic Transform (SHT). Each time-step requires forward and inverse SHTs, and both passes depend on global pencil transposition that redistribute multi-dimensional arrays across compute nodes. At large node counts these global collectives dominate wall-clock time. We investigate GPU-resident lossy compression, using representative fields from the DYAMOND highresolution operational dataset as input, and combining measured GPU compression throughput with SimGrid network simulation, we show that ZFP at 16 bits per value (rate-16), about the same storage budget as float16, matches the communication-time reduction of float16 truncation while delivering approximately 1600× lower mean relative error. ZFP at 8 bits per value (rate8) achieves approximately 1.93× the speedup of float16 while retaining 4× lower mean relative error. Index Terms—spherical harmonic transform, lossy compression, GPU computing, numerical weather prediction, pseudospectral models, all-to-all communication, mixed precision, highperformance computing
I. I NTRODUCTION Numerical weather prediction (NWP) and global climate modeling perform some of the largest and most time-critical simulations in computational science. One of the dominant operational approach to global forecasting is the pseudospectral discretization, exemplified by the ECMWF Integrated Forecasting System (IFS), which today remains a CPU-only production system. GPU ports of IFS and of the wider ECMWF modeling infrastructure are under active development, but as of the most recent ECMWF report they are not yet used operationally [1]. We focus on pseudospectral models specifically because they are communication-bound relative to the grid-point (nearestneighbor) dynamical cores used by other models, such as ICON [2]. The Spherical Harmonic Transform (SHT) is at the heart of a pseudospectral model and it requires two global all-toall-like data redistributions per time-step (Section II-B). This makes inter-node communication, rather than local compute, This work was partially funded under the NRRP, Mission 4 Component 2 Investment 1.4, by the European Union – NextGenerationEU (proj. no. CN 00000013, CUP no. E63C22000970007).
SC26 Workshops, November 15-20, 2026, Chicago, Illinois, USA 979-8-3195-1221-5/26/$31.00 ©2026 IEEE
the eventual bottleneck to GPU-accelerated pseudospectral forecasting and it is precisely this bottleneck that motivates the present study. This pseudospectral core, common to IFS and related models, employs the Spherical Harmonic Transform (SHT) to compute spatial derivatives and to move between spectral coefficients and physical-space grid values. Because the SHT is based on a Gauss–Legendre quadrature in latitude and a Fourier decomposition in longitude, each time-step requires a forward and an inverse SHT. On a distributedmemory machine these transforms involve two global data redistribution, so-called pencil transpositions, that rearrange multi-dimensional arrays via all-to-all collectives: one by rows (from z-pencils to λ-pencils before the longitude FFT) and one by columns (from λ-pencils to θ-pencils before the Legendre transform)(see Section II-B for details). These pencil transpositions scale poorly on modern clusters. While GPU floating-point throughput has grown by orders of magnitude over the past decade, inter-node network bandwidth and latency have not kept pace [3], [4]. As the number of nodes grows, the time spent in global all-to-all communication becomes the dominant contributor to wall-clock time, limiting strong-scaling efficiency and eroding the benefit of faster local computation. A straightforward mitigation is to halve message sizes by casting exchanged buffers to half-precision float (float16) before transmission and restoring float32 on the receiver. This is essentially a form of mixed-precision communication and has attracted growing interest in the NWP community [5]. Its appeal is simplicity: no external library, negligible cast overhead, and a guaranteed 2× reduction in message size. However, the resulting truncation error is uncontrolled, cannot be bounded without knowledge of the data, and may differ dramatically across atmospheric variables and vertical levels. Float16 truncation can introduce errors that accumulate across time-steps and degrade forecast accuracy. We propose online GPU-resident lossy compression as a principled alternative. In our approach, each MPI rank compresses its outgoing communication buffer on the GPU immediately before the all-to-all, transmits the compressed byte stream (of fixed, predictable size), and decompresses on the receiving side before any local computation. By performing compression
and decompression on the GPU, we avoid costly host–device memory transfers. By using ZFP in fixed-rate mode, we obtain deterministic output sizes compatible with standard MPI semantics while exploiting spatial correlations in the atmospheric data to deliver substantially tighter error bounds than float16 for the same storage budget. This paper makes the following contributions: 1) Characterization of NWP data compressibility: We Figure 1: Pencil decomposition of a 3-D data array across MPI analyze representative DYAMOND [6] float32 fields and ranks. Each rank owns a contiguous slab (pencil) along one axis. A characterize the relative error achieved by ZFP fixed-rate pencil transposition redistributes ownership among ranks, realigning compression at multiple rates, comparing against float16 the pencils with a different axis so that the next transform (FFT or Legendre) can proceed locally. The redistribution is implemented truncation. 2) Feasibility study of compressed pencil communication: as an all-to-all collective by rows/columns in which every rank simultaneously sends to and receives from every other rank. To our knowledge, this is the first study to estimate the benefit of GPU-resident compressed communication for a pseudospectral climate/NWP code, which today into a longitude transform (a standard FFT) and a latitude transruns on CPU only. We implement a realistic benchmark form based on the associated Legendre polynomials Pℓm (cos θ). of the SHT’s first all-to-all transposition, measuring The Gauss–Legendre quadrature scheme we consider uses GPU compression and decompression throughput on real Gauss–Legendre latitudes to achieve exact quadrature for hardware and modeling end-to-end communication time polynomials up to the truncation wavenumber; this is the with SimGrid [7] on a Dragonfly topology at node counts scheme used in operational models such as the ECMWF from 4 to 121, representative of current operational NWP IFS, as distinct from equispaced schemes that are more deployments. common in data-analysis applications. In practice, atmospheric 3) Quantitative comparison with float16: We show that state variables are held on a three-dimensional grid with Nλ ZFP rate-16 matches the communication-time reduction longitude points, Nθ Gauss–Legendre latitude points, and Nz of float16 while delivering a mean relative error of 2.5 × vertical levels. 10−7 versus 4×10−4 for float16 (a 1600× improvement) and that ZFP rate-8 achieves 1.93× the speedup of float16 B. Pencil Transpositions On a distributed-memory system, the 3-D data array is at a mean relative error of 10−4 (a 4× improvement). decomposed across MPI ranks. In the pencil decomposition, 4) Real-hardware validation: We corroborate the simulated the domain is partitioned along two of its three dimensions, communication-time trend with measurements on the so that each rank holds a pencil aligned with the third, as CPU partition of the CRESCO8 cluster, grounding the illustrated in Figure 1. Because the FFT (in longitude) and the study in a currently-deployed, production-representative Legendre transform (in latitude) must each operate along their system (Sections IV-E and V-B). respective dimension without distributed data, the transform The remainder of the paper is organized as follows. Section requires two global redistributions: II provides background on the Spherical Harmonic Trans1) All-to-all by rows (transpose from z-pencils to λform, pencil transpositions, and ZFP compression. Section pencils): each rank sends one horizontal slice to every III reviews related work on SHT libraries, distributed FFTs, other rank, enabling the longitude FFT. mixed precision in NWP, and scientific data compression. 2) All-to-all by columns (transpose from λ-pencils to θSection IV describes our experimental methodology, including pencils): a second collective enables the latitude Legendre data selection, benchmark design, and performance modeling. transform. Section V presents our main findings on compressibility and communication time reduction. Some future directions are Both transpositions are local (row-wise or column-wise) all-todiscussed in Section VI. Finally, Section VII summarizes our all calls in which every rank is both a sender and a receiver. In this paper we study the first transposition (rows), noting that the contributions. second is structurally analogous. The communication bottleneck II. BACKGROUND arises because collective bandwidth scales with the number of communicating ranks while per-link network bandwidth does A. Spherical Harmonic Transform not. On clusters where GPU compute throughput is plentiful The Spherical Harmonic Transform decomposes a function but inter-node bandwidth is limited, a gap that continues to defined on the sphere into a linear combination of spherical widen, the all-to-all dominates wall time at moderate-to-large harmonics Yℓm (θ, λ), where θ is the co-latitude, λ is the node counts [3], [4]. longitude, and ℓ, m are the degree and order. For numerical weather prediction, the SHT provides a spectrally accurate way C. ZFP Fixed-Rate Compression to evaluate horizontal derivatives and to apply spectral filters or ZFP [9], [10] is an open-source compressor designed for semi-implicit solvers [8] Computationally, the SHT decomposes multi-dimensional floating-point arrays that exhibit spatial
correlation, as is typical of regularly sampled continuous fields. GPU and 2DECOMP&FFT uses CPU) and heFFTe (GPUZFP partitions the input array into non-overlapping 4d blocks enabled, but MFFT and heFFTe have a different data layout (where d is the array dimensionality) and compresses each and heFFTe needs three communication steps as opposed to block independently by applying a reversible lifting transform two), using compression to reduce message sizes. However, the followed by a precision reduction. ZFP offers three compression comparison is somewhat unbalanced as MFFT uses GPUs for modes. The accuracy mode guarantees that every reconstructed the computation while 2DECOMP&FFT and heFFTe use CPUs. value is within a user-specified absolute error tolerance, but As for heFFTe, the comparison is not entirely fair as heFFTe produces variable-size output that is incompatible with standard uses a different data layout (blocks as opposed to pencils) MPI pre-allocation and is not well-suited to GPU parallelism. and needs three communication steps as opposed to two, and The precision mode targets a fixed number of uncompressed bits communication is known to be the bottleneck for heFFTe per value. The fixed-rate mode assigns exactly r compressed [21]. The key difference between this work and MFFT is that bits per value, producing output of deterministic size and MFFT focuses on FFT and uses a custom lossy compression enabling trivially parallel compression of independent blocks scheme, while we focus on SHT and uses ZFP. Our work is on a GPU. For our initial investigation, we chose to use ZFP in preliminary and a comparison with MFFT is left for future fixed-rate mode. Other compression libraries for scientific data work. Cayrols et al. [24] have studied the speedup of float16 offer different tradeoffs in their accuracy models, performance truncation during communication for FFT, while they achieve profiles, and GPU support (see Section III-E). Those will be a significant speedup, the approach is naive and not tunable, explored in a future work. as the error cannot be controlled. Furthermore, the accuracy gain of using float32 computation with float16 communication III. R ELATED W ORK over using both computation and communication in float16 may not justify this approach. A. SHT Libraries for Distributed Systems Several open-source SHT libraries are available, but few C. Mixed Precision in NWP support both MPI distribution and GPU acceleration simultaThe use of reduced floating-point precision in atmospheric neously. SHTns [11] and DUCC [12] are highly optimized for models has received considerable attention. The ECMWF IFS shared-memory machines (CPU+GPU) but do not support MPI. now runs parts of its computation in single precision [5]. Our HEALPix [13], CHarm [14] and related astronomical codes work is complementary: rather than reducing the precision of operate on a different pixelization scheme and typically analyze arithmetic, we reduce the precision of communication buffers, a single scalar field, making them unsuitable for NWP with using a compressor that adapts its bit-allocation to local field hundreds of vertical levels. ISPACK [15] and Libsharp [16] structure. are consolidated libraries that have MPI support but lack GPU support. The ECMWF ecTrans library [17] is, to our knowledge, D. Compression of NWP Data the only library suitable for large-scale operational NWP: it Various studies have investigated the use of lossy comsupports MPI distribution, operates on Gauss–Legendre grids, pression for NWP data, but they have focused on archival and has been partially ported to GPU. However, as of the time compression of large datasets, not on online compression of of writing, operational ECMWF forecasts are still produced communication buffers. Tintó Prims et al. [25] study the effect on CPUs, and the GPU port has not yet been fully optimized of lossy compression of whole datasets for analysis purposes for production use [1] with a custom compressor. Huang et al. [26] study the feasibility of custom ML compressors for archival purposes. Those works B. Distributed FFT and Pencil Communication exploit temporal correlations across many time-steps and long The pencil transposition pattern is essentially identical to records and are unsuitable for online per-message compression. that used in distributed three-dimensional FFTs. Libraries such Our setting is fundamentally different: we compress a small as 2DECOMP&FFT [18], [19], P3DFFT [20], and heFFTe spatial patch of a single time-step, in real time, on a GPU, [21] implement pencil transpositions for FFTs, and their with no access to past or future fields. Klöwer et al. [27] do a communication structure and scaling behavior are directly study on the effect of lossy compression on the CAMS dataset relevant to the SHT setting. None of these libraries currently (float64) using ZFP as the present work. However, their work support variable-length compressed messages or expose a focuses on the effect of lossy compression on accuracy and not compression hook in their communication pipeline. 2DE- on the communication time, which is the focus of our work. COMP&FFT [18], [19] exposes an API for general pencil They study the compression of a single time-step, without transposition and cuDecomp [22] is a library focused on GPU- knowledge of past and future data, but they still compress accelerated pencil transpositions with an automatic choice data for the whole 3D atmospheric fields. In our work we of communication strategy and backend. A closely related are bound to compression on data resident on a single GPU, contribution to our work is the MFFT library [23], which and so the achievable compression ratio is limited. Their work applies lossy compression to the inter-process communication however is complementary to ours, and we plan to extend our buffers of a distributed FFT. MFFT reports performance work by using their analysis of the effect of lossy compression improvements over 2DECOMP&FFT (however MFFT uses on the forecast accuracy to inform the choice of compression
parameters for the best tradeoff between communication time and forecast accuracy.
B. Compression Configuration
We use ZFP version 1.0.1 in fixed-rate mode. Fixed-rate mode was selected for an initial performance investigation E. Compression of Scientific Floating-Point Data because it is the only mode with GPU support and because it Many compressors exist for scientific data, including ZFP [9], produces messages of deterministic size, which is compatible SZ and its variants [28], [29], and MGARD [30]. See also Di with standard MPI semantics and easier to integrate into et al. [31] for a survey. Compressor in their accuracy models, existing SHT libraries. Future works may explore accuracyperformance profiles, and GPU support. For our feasibility mode compression, which offers tighter error control at the study the choice of compressor is secondary; ZFP was selected cost of variable message sizes. We evaluate two rates. Rate-16 because its CUDA backend (cuZFP) is mature, its fixed- (16 bits per value) produces messages of the same size as rate mode produces outputs of predictable size, it has been float16 truncation, enabling a direct comparison of accuracy extensively validated in HPC settings [10], and has already at equal communication cost. Rate-8 (8 bits per value) halves been used for compression of atmospheric data [27]. The results message size again, at the cost of higher compression error. we report are therefore a conservative lower bound on what a compressor better adapted to atmospheric field structure could C. GPU Compression Throughput achieve. Compression and decompression were performed on a NVIDIA A100 SXM4 GPU (19.5 TFlop/s of theoretical F. Compressed Communication Middleware float32 performance) with 40 GB of memory, using the cuZFP Recent works on compression-aware MPI and network CUDA backend. We measured wall-clock throughput by using middleware [32]–[35] apply generic compression algorithms at the cuZFP CUDA backend. We measured wall-clock throughput the communication layer. Such approaches do not necessitate by compressing and decompressing the message buffers for to modify application code, but they are not suitable for our use representative fields T and O . Measurements were taken after 3 case. These compressors are typically data-agnostic and operate a warm-up phase to exclude GPU initialization overhead. on flattened one-dimensional byte buffers, which prevents them CPU compression throughput was also measured on a AMD from exploiting the multi-dimensional block structure on which EPYC 7742 64-Core Processor ( 2.5 GHz of clock frequency) accuracy guarantees of scientific compression libraries depend. for comparison. Compression is compute-bound and multi-core Because MPI primitives do not expose multi-dimensional array CPU does not provide the necessary parallelism for efficient metadata to the communication layer, injecting a structure- compression. We present this result quantitatively in Section aware compressor requires either application-level wrapping V-D to motivate the GPU-only approach. (our approach) or substantial MPI implementation changes. D. Network Simulation with SimGrid IV. M ETHODOLOGY We use state-of-the-art SimGrid simulator [7] with its SMPI A. Experimental Data module to model all-to-all communication time at scale. Its key We use output from the DYAMOND [6] high-resolution advantage for our purposes is that it eliminates measurement global storm-resolving simulation at 2560×1281 grid points, ap- noise arising from production network congestion, enabling proximately twice the horizontal resolution of current ECMWF clean comparisons between communication strategies across a operational forecasts. The dataset is distributed in 32-bit wide range of node counts without requiring access to a large single-precision floating point, making it representative of the cluster. We modeled a Dragonfly topology with the following precision used in modern NWP systems (operational models are parameters: 27 groups with 4 group-group links, 16 chassis increasingly moving from float64 to float32 [5]). Operational per group with 3 chassis-to-chassis links, 4 routers per chassis NWP fields fall into two distinct categories. Dynamical- with 1 router-to-router link, and 1 nodes per router, yielding core variables such as temperature and ozone are stored and 6480 nodes in total. Link bandwidth was set to 200 Gbps communicated as native float32. Many tracer and diagnostic with a latency of 1 µs per hop. This topology is similar in variables, however, are stored with an integer quantization scale and configuration to the JUPITER system. Node counts scheme (an integer array paired with a scale factor and offset), ranged from 4 to 121 nodes (as listed in Table I), a range which reduces their dynamic range and makes them unsuitable directly relevant to current operational NWP deployments. We test cases for a floating-point compressor: the piecewise- chose this range rather than a larger future-scale simulation constant quantization creates sharp step discontinuities that to demonstrate that the compression approach is beneficial for are pathological for block-transform compressors such as present-day routine weather forecasts, not only at hypothetical ZFP. We therefore focus on temperature (T ) and ozone (O3 ), exascale node counts. which are representative of the float32 dynamical-core variables The simulation proceeds as follows. For each node count that dominate communication volume and forecast sensitivity. N and each compression strategy, SimGrid simulates the We focus on compressing floating point variables as we are communication collective using the message sizes from Table I. interested in the performance of compressing data for numerical The compression and decompression overheads, are added as prediction (online), as opposed to compressing data for archival a compute-side delay on the sender (compression) and receiver purposes (offline). (decompression) before and after the collective, respectively.
Float16 truncation is modeled with negligible cast overhead (the cast from float32 to float16 is a GPU register operation). Total time is reported as tcompress +tcommunication +tdecompress , where tcommunication is the simulated communication time and the remaining terms are measured. This serial model is a conservative upper bound, as compression and communication could in principle be pipelined per block. E. CRESCO8 Real-Hardware Measurements To complement the SimGrid simulation, we measured communication time directly on the CPU partition of the CRESCO8 cluster (ENEA, Portici, Italy). CRESCO8 CPU partition comprises 760 compute nodes (dual-socket, 64-core Intel Xeon Platinum 8592+, 1.9–3.8 GHz, 512 GB RAM) Nodes are connected through a 200 Gb/s Mellanox NDR InfiniBand fabric in a Dragonfly topology, the same topology class used in the SimGrid model. We ran the row-wise all-to-all transposition on CRESCO8’s CPU nodes and measured communication time for the four strategies (float32, float16, ZFP rate-16, ZFP rate-8) at world sizes of 4, 16, and 64. Operational NWP codes run on CPU only, so compression and decompression segments are not measured natively; they reuse the NVIDIA A100 throughput of Section IV-C, while only the communication segments were measured on CRESCO8.
V. R ESULTS A. SimGrid Simulation of Communication Time Figure 2 shows stacked bar plots of total time (compression overhead + communication time + decompression overhead) as a function of node count N , for four strategies: baseline float32, float16 truncation, ZFP rate-16, and ZFP rate-8. The following observations summarize the results. Rate-16 vs. float16. At small node counts (N ≤ 25), GPU compression overhead is non-negligible relative to the communication time, and rate-16 is slightly slower than float16. At N ≈ 36–100 nodes, the compression overhead becomes small relative to the reduced communication time and rate-16 matches float16. At large node counts (N = 121), the message size of float16 and rate-16 is comparable and compression overhead becomes negligible, so the speedup of both strategies is similar. This matched speedup comes at a much lower error: as Section V-C shows, rate-16’s mean relative error is 1600× lower than float16’s (Figure 4). Rate-8. Rate-8 produces messages about half the size of rate-16, and at large node counts (N ≥ 100) it achieves approximately 1.93× the speedup of float16 over the float32 baseline. At all node counts, rate-8 is faster than float16. As Section V-C shows, rate-8 also retains a 4× lower mean relative error than float16 (Figure 5).
F. Error Metric As a proxy for the error on the downstream application, we report the relative error: ϵ=
∥x − x̂∥ , ∥x∥
where x are the original values and x̂ are the reconstructed values. For float16 truncation the analogous maximum relative error is computed by casting the float32 values to float16 and back and applying the same formula.
N (nodes)
Float32 (MB)
Float16 (MB)
Rate-16 (MB)
Rate-8 (MB)
4 9 16 25 36 49 64 81 100 121
212.50 62.38 26.56 13.50 7.62 4.81 3.32 2.31 1.62 1.23
106.25 31.19 13.28 6.75 3.81 2.41 1.66 1.15 0.81 0.62
106.25 33.54 14.06 7.00 4.23 2.58 1.95 1.25 1.00 0.62
53.13 16.77 7.03 3.50 2.12 1.29 0.98 0.62 0.50 0.31
Table I: MPI message sizes per rank-pair (in MB) for the all-to-all-byrows transposition as a function of node count N . Sizes are computed for the DYAMOND 2560 × 1281 grid. Float32 is the uncompressed baseline; float16 and ZFP rate-16 are similar in size by construction (both store 16 bits per value); ZFP rate-8 approximately halves that again.
Figure 2: Simulated end-to-end time (compression + communication + decompression) for the four strategies, as a function of node count N from 4 to 121, using GPU compression and decompression throughput measured on the NVIDIA A100 (Section IV-C). The compression and decompression overhead is a small and shrinking fraction of the total as N grows, becoming almost negligible at the largest node counts.
B. CRESCO8 Communication Measurements Figure 3 shows stacked bar plots of total time (compression overhead + communication time + decompression overhead) on the CRESCO8 CPU partition (Section IV-E) for the four strategies, as a function of MPI world size. As in the simulated results of Figure 2, the compression/decompression overhead shrinks relative to communication as scale grows. The real measurements reach world size 64. The following observations summarize the results. Rate-16 vs. float16. At N = 4, compression overhead dominates the total time, though this world size is too small to be representative of realistic
deployments. At N = 16, the compressed message is faster than the uncompressed baseline but not yet competitive with float16. At N = 64, rate-16 is slightly slower than float16; however, as Section V-C shows, rate-16’s mean relative error is 1600× lower than float16’s (Figure 4). Rate-8. At N = 4, compression overhead again dominates the total time, and this world size is not representative of realistic deployments. At N = 16, rate-8 is faster than float16 while retaining better accuracy. At N = 64, rate-8 is considerably faster than float16, and, as Section V-C shows, retains a 4× lower mean relative error (Figure 5).
Figure 4: Relative error distributions for ZFP rate-16 and float16 truncation applied to the DYAMOND summer temperature field (201608-01 00Z). ZFP rate-16 achieves a mean relative error of 2.5 × 10−7 ; float16 truncation achieves a mean of 4×10−4 . Compression achieves a 1600× lower mean relative error than float16 at the same storage cost. Ozone (O3 ) shows a very similar relative-error distribution and is omitted for brevity.
Figure 3: Measured communication time on the CRESCO8 CPU partition for the four strategies, as a function of MPI world size (4, 16, 64). Compression and decompression segments overlay the NVIDIA A100 throughput of Section IV-C; only the communication segments were measured on CRESCO8 itself.
C. Data Compressibility and Compression Error Table I shows how individual message sizes decrease as node count grows. Message sizes have been calculated on the DYAMOND data for the temperature field. Figures 4 and 5 show the distributions of relative error for temperature (T ) over all 43 ZFP blocks in the DYAMOND field, comparing Figure 5: Relative error distributions for ZFP rate-8 and float16 truncation applied to the same DYAMOND temperature field. ZFP ZFP compression against float16 truncation at rate-16 and rate- rate-8, which uses around 8 bits per value (half the storage of float16), 8 respectively. Figure 4 shows a clear difference at similar achieves a mean relative error of 10−4 , better than the float16 mean storage cost. ZFP rate-16 achieves a mean relative error of of 4 × 10−4 . As with rate-16, ozone (O3 ) shows a very similar 2.5 × 10−7 , compared to 4 × 10−4 for float16 truncation, a relative-error distribution and is omitted for brevity. factor of 1600× lower. Figure 5 shows that ZFP rate-8, using half the storage of float16, achieves a mean relative error of 10−4 , within a factor of 4 of float16’s mean relative error, while approximately halving the message size again. This (compression overhead + communication time + decompression error analysis suggests that for a typical atmospheric field, overhead) as a function of node counts N . At low-to-moderate ZFP’s block-transform strategy achieves an acceptable error node counts (N ≤ 36), the data volume per node is large and with a better accuracy–performance trade-off than float16. This CPU compression time exceeds by one order of magnitude the may not be the case for extreme events, but this has to be communication savings from reduced message size, making the further investigated in future work. We are aware that the error approach counter-productive. At large node counts (N = 121), analysis is limited to a couple of atmospheric variables and a message sizes are small enough that CPU compression becomes single snapshot. We plan to extend our work by using different feasible in principle, but large-scale NWP is a weak-scaling fields from various weather condition following the analysis workload: as node counts increase, the problem size typically increases proportionally, so the data volume per node remains of Klöwer et al. [27]. approximately constant. This means that CPU compression D. CPU Compression: Why It Is Not Viable time does not decrease as N grows in a realistic production Figure 6 shows a stacked bar chart similar to the one in scenario, confirming that GPU compression is the only viable Figure 2 with the timing breakdown for CPU compression pathway.
(see Section III-E), and custom-designed compressors tailored to atmospheric fields may achieve substantially better accuracy at the same rate. Our results are therefore a conservative lower bound on the achievable accuracy improvement. We reported error at the communication buffer level. The effect of this error on the final spherical harmonic coefficients and, downstream, on forecast accuracy over multi-day integrations requires a full end-to-end experiment, see Klöwer et al. [27] for a preliminary analysis of the effect of lossy compression on forecast accuracy. VII. C ONCLUSION
Figure 6: Same breakdown as Figure 2, using CPU compression and decompression throughput (AMD EPYC 7742) instead of GPU throughput. CPU compression time exceeds the communication time for most node counts, only becoming comparable to it at the largest node counts.
E. Accuracy–Performance Trade-off Table II summarizes the accuracy–performance trade-off at N = 121 nodes, representative of a large-scale operational run. Float16 truncation achieves a 2.05× speedup over float32 at a mean relative error of 4 × 10−4 . ZFP rate-16 achieves roughly the same speedup (1.95×) at a mean relative error of 2.5 × 10−7 , a factor of 1600× lower at a negligible additional cost. ZFP rate-8 achieves a 4.1× speedup at a mean relative error of 10−4 , which is 4× better than float16 at nearly double the speedup. These results establish a clear improvement over float16 truncation. ZFP rate-16 achieves the same communication-time reduction at a negligible additional cost, while delivering a mean relative error that is three orders of magnitude lower. ZFP rate-8 achieves an even better speedup at a mean relative error that is still better than float16 error. VI. F UTURE D IRECTIONS In this work we only addressed the first all-to-all transposition (by rows). The second transposition (by columns) is structurally analogous; the compression rationale and error analysis apply equally. End-to-end SHT performance would require both transpositions to be compressed, and the combined speedup would be somewhat smaller than the per-transposition figures reported here, as local compute (FFT, Legendre transform) is not accelerated. ZFP is a strong but general-purpose compressor. Other floating point compression libraries exist Table II: Summary of the accuracy–performance trade-off at N = 121 nodes. Mean relative error is reported over all 43 ZFP blocks of the temperature field. Speedup is relative to the float32 baseline. Best values are underlined. Strategy Float32 (baseline) Float16 ZFP rate-16 ZFP rate-8
Speedup ↑
Mean rel. error ↓
1.0× 2.05× 1.95× 3.97×
— 4 × 10−4 2.5 × 10−7 1 × 10−4
We have presented a study on GPU-resident lossy compression as an alternative to float16 truncation for reducing the cost of MPI all-to-all pencil transpositions in the Spherical Harmonic Transform. Using measured GPU compression throughput from the cuZFP backend and a SimGrid network simulation on a Dragonfly topology at node counts representative of current operational NWP deployments, we have shown that ZFP fixed-rate compression at 16 bits per value achieves roughly the same communication-time reduction as float16 truncation while delivering a mean relative error of 2.5 × 10−7 , a factor of 1600× lower than float16’s 4 × 10−4 on the DYAMOND temperature field. ZFP at 8 bits per value achieves approximately 1.93× the speedup of float16 at a mean relative error of 10−4 , which is 4× better than float16 at nearly double the speedup. We corroborated the trend of these simulated results with real measurements on the CPU partition of the CRESCO8 cluster, grounding the study in a currentlydeployed, production-representative system. Works like [25], [27] suggest that achievable compression ratio for atmospheric data may be higher than what we claimed, and so the accuracy– performance trade-off may be even better than what we report here. However, this has to be further investigated in future work balancing compression ratio, compression time and downstream forecast accuracy. These results establish that the accuracy– performance trade-off curve of ZFP fixed-rate strictly dominates that of float16 truncation. We have also shown that CPU-based compression is not viable as a drop-in replacement, confirming that GPU-resident compression is essential. Operational pseudospectral climate and NWP codes run on CPUs today, and GPU-resident compressed communication is, to our knowledge, a currently-unexplored option for them. As these codes migrate to GPU, we argue this is an opportunity the climate/NWP community should weigh alongside the engineering cost of the migration itself. ACKNOWLEDGEMENTS This work was partially funded under the NRRP, Mission 4 Component 2 Investment 1.4, by the European Union – NextGenerationEU (proj. no. CN 00000013, CUP no. E63C22000970007). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or The European Research Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.
R EFERENCES [1] ECMWF, “ECMWF Newsletter 186 - Winter 2025/26,” 2026. [2] D. Klocke, C. Frauen, J. F. Engels, D. Alexeev, R. Redler, R. Schnur, H. Haak, L. Kornblueh, N. Brüggemann, F. Chegini, M. Römmer, L. Hoffmann, S. Griessbach, M. Bode, J. Coles, M. Gila, W. Sawyer, A. Calotoiu, Y. Budanaz, P. Mazumder, M. Copik, B. Weber, A. Herten, H. Bockelmann, T. Hoefler, C. Hohenegger, and B. Stevens, “Computing the Full Earth System at 1km Resolution,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’25. New York, NY, USA: Association for Computing Machinery, Nov. 2025, pp. 125–136. [3] D. De Sensi, L. Pichetti, F. Vella, T. De Matteis, Z. Ren, L. Fusco, M. Turisini, D. Cesarini, K. Lust, A. Trivedi, D. Roweth, F. Spiga, S. Di Girolamo, and T. Hoefler, “Exploring GPU-to-GPU Communication: 2024 International Conference for High Performance Computing, Networking, Storage and Analysis, SC 2024,” Sc24: International Conference for High Performance Computing, Networking, Storage and Analysis, pp. 1–15, 2024. [4] D. Unat, I. Turimbetov, M. K. T. Issa, D. Sağbili, F. Vella, D. De Sensi, and I. Ismayilov, “The Landscape of GPU-Centric Communication,” 2024. [5] F. Váňa, P. Düben, S. Lang, T. Palmer, M. Leutbecher, D. Salmond, and G. Carver, “Single Precision in Weather Forecasting Models: An Evaluation with the IFS,” Monthly Weather Review, vol. 145, no. 2, pp. 495–502, Feb. 2017. [6] B. Stevens, M. Satoh, L. Auger, J. Biercamp, C. S. Bretherton, X. Chen, P. Düben, F. Judt, M. Khairoutdinov, D. Klocke, C. Kodama, L. Kornblueh, S.-J. Lin, P. Neumann, W. M. Putman, N. Röber, R. Shibuya, B. Vanniere, P. L. Vidale, N. Wedi, and L. Zhou, “DYAMOND: The DYnamics of the Atmospheric general circulation Modeled On Nonhydrostatic Domains,” Progress in Earth and Planetary Science, vol. 6, no. 1, p. 61, Sep. 2019. [7] H. Casanova, A. Giersch, A. Legrand, M. Quinson, and F. Suter, “Lowering entry barriers to developing custom simulators of distributed applications and platforms with SimGrid,” Parallel Computing, vol. 123, p. 103125, Mar. 2025. [8] J. Coiffier, Fundamentals of Numerical Weather Prediction, 1st ed. Cambridge University Press, Dec. 2011. [9] P. Lindstrom, “Fixed-Rate Compressed Floating-Point Arrays,” IEEE Transactions on Visualization and Computer Graphics, vol. 20, no. 12, pp. 2674–2683, Dec. 2014. [10] J. Diffenderfer, A. L. Fox, J. A. Hittinger, G. Sanders, and P. G. Lindstrom, “Error Analysis of ZFP Compression for Floating-Point Data,” SIAM Journal on Scientific Computing, vol. 41, no. 3, pp. A1867–A1898, Jan. 2019. [11] N. Schaeffer, “Efficient Spherical Harmonic Transforms aimed at pseudospectral numerical simulations,” https://arxiv.org/abs/1202.6522v5, Feb. 2012. [12] M. Reinecke, “DUCC: Distinctly Useful Code Collection,” Astrophysics Source Code Library, p. ascl:2008.023, Aug. 2020. [13] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, “HEALPix: A Framework for HighResolution Discretization and Fast Analysis of Data Distributed on the Sphere,” The Astrophysical Journal, vol. 622, no. 2, p. 759, Apr. 2005. [14] B. Bucha, “CHarm: C/Python library for spherical harmonic transforms in planetary geodesy,” Earth Science Informatics, vol. 19, no. 3, p. 29, Mar. 2026. [15] K. Ishioka, “A New Recurrence Formula for Efficient Computation of Spherical Harmonic Transform,” Journal of the Meteorological Society of Japan. Ser. II, vol. 96, no. 2, pp. 241–249, 2018. [16] M. Reinecke and D. S. Seljebotn, “Libsharp – spherical harmonic transforms revisited,” Astronomy & Astrophysics, vol. 554, p. A112, Jun. 2013. [17] ECMWF, “IFS Documentation CY49R1 - Part VI: Technical and Computational Procedures,” 2024. [18] N. Li and S. Laizet, “2decomp & fft-a highly scalable 2d decomposition library and fft interface,” in Cray User Group 2010 Conference, 2010, pp. 1–13. [19] S. Rolfo, C. Flageul, P. Bartholomew, F. Spiga, and S. Laizet, “The 2DECOMP&FFT library: An update with new CPU/GPU capabilities,” Journal of Open Source Software, vol. 8, no. 91, p. 5813, Nov. 2023.
[20] D. Pekurovsky, “P3DFFT: A framework for parallel computations of Fourier transforms in three dimensions,” SIAM Journal on Scientific Computing, vol. 34, no. 4, pp. C192–C209, Jan. 2012. [21] A. Ayala, S. Tomov, A. Haidar, and J. Dongarra, “heFFTe: Highly Efficient FFT for Exascale,” in Computational Science – ICCS 2020, V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos, and J. Teixeira, Eds. Cham: Springer International Publishing, 2020, pp. 262–275. [22] J. Romero, P. Costa, and M. Fatica, “Distributed-memory simulations of turbulent flows on modern GPU systems using an adaptive pencil decomposition library,” in Proceedings of the Platform for Advanced Scientific Computing Conference, ser. PASC ’22. New York, NY, USA: Association for Computing Machinery, Jul. 2022, pp. 1–11. [23] Y. Zhao, F. Liu, W. Ma, H. Li, Y. Peng, and C. Wang, “MFFT: A GPU Accelerated Highly Efficient Mixed-Precision Large-Scale FFT Framework,” ACM Transactions on Architecture and Code Optimization, vol. 20, no. 3, pp. 1–23, Sep. 2023. [24] S. Cayrols, J. Li, G. Bosilca, S. Tomov, A. Ayala, and J. Dongarra, “Lossy all-to-all exchange for accelerating parallel 3-D FFTs on hybrid architectures with GPUs,” in 2022 IEEE International Conference on Cluster Computing (CLUSTER). Heidelberg, Germany: IEEE, Sep. 2022, pp. 152–160. [25] O. Tintó Prims, R. Redl, M. Rautenhaus, T. Selz, T. Matsunobu, K. R. Modali, and G. Craig, “The effect of lossy compression of numerical weather prediction data on data analysis: A case study using enstoolscompression 2023.11,” Geoscientific Model Development, vol. 17, no. 24, pp. 8909–8925, Dec. 2024. [26] L. Huang and T. Hoefler, “Compressing multidimensional weather and climate data into neural networks,” Apr. 2023. [27] M. Klöwer, M. Razinger, J. J. Dominguez, P. D. Düben, and T. N. Palmer, “Compressing atmospheric data into its real information content,” Nature Computational Science, vol. 1, no. 11, pp. 713–724, Nov. 2021. [28] X. Liang, K. Zhao, S. Di, S. Li, R. Underwood, A. M. Gok, J. Tian, J. Deng, J. C. Calhoun, D. Tao, Z. Chen, and F. Cappello, “SZ3: A Modular Framework for Composing Prediction-Based Error-Bounded Lossy Compressors,” IEEE Transactions on Big Data, vol. 9, no. 2, pp. 485–498, Apr. 2023. [29] J. Tian, S. Di, K. Zhao, C. Rivera, M. H. Fulp, R. Underwood, S. Jin, X. Liang, J. Calhoun, D. Tao, and F. Cappello, “cuSZ: An Efficient GPU-Based Error-Bounded Lossy Compression Framework for Scientific Data,” in Proceedings of the ACM International Conference on Parallel Architectures and Compilation Techniques, ser. PACT ’20. New York, NY, USA: Association for Computing Machinery, Sep. 2020, pp. 3–15. [30] Q. Gong, J. Chen, B. Whitney, X. Liang, V. Reshniak, T. Banerjee, J. Lee, A. Rangarajan, L. Wan, N. Vidal, Q. Liu, A. Gainaru, N. Podhorszki, R. Archibald, S. Ranka, and S. Klasky, “MGARD: A multigrid framework for high-performance, error-controlled data compression and refactoring,” SoftwareX, vol. 24, p. 101590, Dec. 2023. [31] S. Di, J. Liu, K. Zhao, X. Liang, R. Underwood, Z. Zhang, M. Shah, Y. Huang, J. Huang, X. Yu, C. Ren, H. Guo, G. Wilkins, D. Tao, J. Tian, S. Jin, Z. Jian, D. Wang, M. H. Rahman, B. Zhang, S. Song, J. Calhoun, G. Li, K. Yoshii, K. Alharthi, and F. Cappello, “A Survey on ErrorBounded Lossy Compression for Scientific Datasets,” ACM Computing Surveys, vol. 57, no. 11, pp. 287:1–287:38, Jun. 2025. [32] Q. Zhou, P. Kousha, Q. Anthony, K. Shafie Khorassani, A. Shafi, H. Subramoni, and D. K. Panda, “Accelerating MPI All-to-All Communication with Online Compression on Modern GPU Clusters,” in High Performance Computing, A.-L. Varbanescu, A. Bhatele, P. Luszczek, and B. Marc, Eds. Cham: Springer International Publishing, 2022, vol. 13289, pp. 3–25. [33] Q. Zhou, C. Chu, N. S. Kumar, P. Kousha, S. M. Ghazimirsaeed, H. Subramoni, and D. K. Panda, “Designing High-Performance MPI Libraries with On-the-fly Compression for Modern GPU Clusters,” in 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). Portland, OR, USA: IEEE, May 2021, pp. 444–453. [34] J. Huang, S. Di, X. Yu, Y. Zhai, Z. Zhang, J. Liu, X. Lu, K. Raffenetti, H. Zhou, K. Zhao, Z. Chen, F. Cappello, Y. Guo, and R. Thakur, “An Optimized Error-controlled MPI Collective Framework Integrated with Lossy Compression,” Jan. 2024. [35] S. Ma, C. L. Lao, Z. Xu, Z. Wang, Z. Mao, D. Meng, J. Zhen, J. Wu, I. Stoica, Y. Wang, and Y. Zhou, “UCCL-Zip: Lossless Compression Supercharged GPU Communication,” Apr. 2026.