Performance report and optimized implementation of Weather & Climate Dwarfs on GPU, MIC and Optalysys Optical Processor

Alan Gray; Alex Macfaden; Cyril Mazauric; David Guibert; Erwan Raffin; Jacob Poulsen; Per Berg; Peter Messmer; Xavier Vigouroux

REVIEW 4 major objections 5 minor 15 references

Performance report and optimized implementation of Weather & Climate Dwarfs on GPU, MIC and Optalysys Optical Processor

T0 review · 4 major / 5 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Restructuring the data flow of weather and climate codes' canonical computation patterns — the 'dwarfs' — yields up to 23x GPU speedups and a first purely optical route to the spherical-harmonics spectral transform.

desk verdict A useful engineering deliverable with a transparent but untransferable GPU baseline; the 23x headline is real but mostly overhead removal, and the optical piece is honestly labeled as work-in-progress. read the letter →

arxiv 1908.06096 v1 pith:KWY3PC6S submitted 2019-08-16 cs.DC

Cyril Mazauric , Erwan Raffin , Xavier Vigouroux , David Guibert , Alex Macfaden , Jacob Poulsen , Per Berg , Alan Gray

show 1 more author

Peter Messmer

This is my paper · ORCID

classification cs.DC

keywords weatherandclimatedwarfsspectraltransformsphericalharmonicsGPUoptimizationrooflinemodelopticalcorrelatorMPDATAcompilerdirectives

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This deliverable report studies the recurring computation-and-communication patterns — called 'dwarfs' — that dominate numerical weather prediction and climate codes, and asks how fast each can run on a single modern node. It claims that compiler directives alone give 1.1x–2.5x on CPUs; that restructuring the GPU code for data locality cuts the spherical-harmonics spectral transform by 23x (to around 40% of the roofline ceiling) and the MPDATA kernel by 9x; and that a purely optical implementation of the spectral transform is feasible, changing the scaling regime from $\mathcal{O}(N^3)$ to $\mathcal{O}(N)$. A reader should care because the spectral transform is the backbone of the spectral formulation of global forecasting models, and the report shows concretely where the headroom sits: after restructuring, most kernels reach 70–90% of achievable memory bandwidth, while FFT library calls and inter-node all-to-all communication remain the binding constraints. The report itself limits all results to single-node problem sizes and defers multi-node and energy comparisons to later work.

What carries the argument

The argument is carried by a small set of named objects. The roofline model is the yardstick: a plot of achievable performance against operational intensity (flops per byte) that classifies each kernel as memory-bound or compute-bound and fixes the ceiling every measured gain is referenced to. The 'dwarf' concept is the unit of generality: each dwarf is a canonical computation-and-communication pattern lifted from operational weather codes, so an optimization proven on the dwarf is meant to transfer to real models. On the GPU side, the mechanism is the restructuring recipe: hoisting allocations out of the timestep loop; collapsing the spherical-dimension and field loops into one collapsed parallel space with masked extents; zero-padding unequal-size matrix multiplications into a single uniform batched call; flattening deep data structures into constant-memory arrays; and tiling transposes so global-memory traffic is coalesced. On the optical side, the mechanism is the optical correlator — a '4f' system (two lenses separated by their focal lengths, with a filter at the shared focal plane) in which a lens performs a physical Fourier transform and the filter multiplies the spectrum, so the whole device evaluates a convolution in the time light takes to cross it; the paper's contribution is a multichannel astigmatic variant whose cylindrical lenses do parallel 1D Fourier transforms along one axis while the correlator stage extracts Legendre coefficients along the other.

What would settle it

Take the same TL159 spectral-transform case on the same GPU and build a control baseline that already hoists allocations out of the timestep loop and batches kernel launches, then compare it with the optimized version: if the mature baseline runs within a small factor of the optimized code, the order-of-magnitude claims are an artifact of the immature baseline, not a property of the algorithm. A second decisive check would compare optically extracted Fourier–Legendre coefficients against a double-precision digital reference at full model resolution, to see whether the optical route can reach the numerical fidelity the report leaves open.

Watch

Extended reading notes

Core claim

The paper's central claim is that the spherical-harmonics spectral transform dwarf's poor GPU performance was never intrinsic to the mathematics. The baseline GPU port ran one to two orders of magnitude below the roofline because its profile was dominated by API calls — per-timestep memory allocations, scattered copies, and launches of a huge number of tiny kernels — and because only two of the three available levels of parallelism (the two spherical dimensions and the field index) were exposed, with a sequential outer loop hiding the third. By hoisting all memory management out of the timestep loop, collapsing the loop nest so all three levels map to the hardware (masking latitude-dependent extents), padding and batching the variable-size matrix multiplications into one uniform batched library call, removing synchronization between FFT calls, replacing deep data-structure accesses with flat arrays that the compiler places in constant memory, and tiling transposes for coalesced access, the same algorithm runs 23x faster on a single GPU and reaches around 40% of the roofline, with matrix multiplications above 80% of peak and most memory-bound kernels at 70–90% of achievable bandwidth; the FFT library calls, at 10–40% of their roofline, are named as the main remaining limiter. The same recipe is reported for the MPDATA flux-divergence kernel, which improved 9.4x to 66% of achievable bandwidth. On the optical side, the paper claims that a coherent optical correlator can extract Fourier–Legendre coefficients as optical dot products, and presents a prototype multichannel astigmatic processor that cascades a 1D optical Fourier-transform stage with a Legendre correlation stage; it demonstrates correlation peaks and derives $\mathcal{O}(N)$ scaling for the transform in place of the digital $\mathcal{O}(N^3)$, while conceding that quantitative coefficient extraction remains challenging and that spatial-light-modulator refresh rates and read-out electronics currently cap throughput.

Load-bearing premise

The load-bearing premise is that the original GPU port — fresh from the compiler, dominated by per-timestep allocations and kernel-launch overhead, and running one to two orders of magnitude below the roofline — is a fair and representative baseline, so the 23x and 9x factors measure the algorithm's headroom rather than the baseline's immaturity.

Editorial extensions

If this is right

If the GPU results hold, the spherical-harmonics spectral transform is not inherently slow on a GPU: after restructuring, most kernels sit at 70–90% of achievable memory bandwidth, and the report names the remaining lever — a custom FFT that batches unequal sizes, which would attack the one component still at 10–40% of its roofline.
The same categories of changes transfer across memory-bound dwarfs: MPDATA's kernel improved 9.4x with loop collapse, flat constant arrays, and replacing division by reciprocal multiplication, so other memory-bound weather and climate kernels are candidates for the same treatment.
A fast spectral-transform coprocessor would not automatically accelerate a full global forecast model: the report's timestep-profile analysis shows that all-to-all communication between nodes dominates the transform phase, leaving post-processing of model output as the natural niche for such a device.
On CPUs the directive-based work is near the limit of what non-intrusive changes buy: gains of 1.1x–2.5x come from vectorization hints, system tuning (huge pages, memory allocator, turbo frequency), and scheduling choices, and the report concludes the dwarfs are memory-bound with further gains requiring deeper data-structure refactoring.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The 23x and 9x figures almost certainly overstate what a mature implementation would gain: a large share of the improvement is removal of API-call and kernel-launch overhead that a well-built GPU port would never have had, so the transferable content is the restructuring recipe, not the factor itself.
A clean test the report leaves implicit: build the proposed batched unequal-size FFT and measure the whole timestep's roofline fraction — the paper's own kernel-level numbers imply the overall 40% could move toward the 70–90% its other kernels already reach.
The optical work is best read as a feasibility and scaling-regime demonstration rather than a performance result: the report's own discussion concedes that multilevel modulator refresh rates near 100 Hz and read-out electronics cap throughput, so the durable claim is architectural — a lens system evaluates the transform's inner-product structure at $\mathcal{O}(1)$ per projection — and the testabl
The report's method of folding a non-uniform sampling grid into the filter rather than interpolating the data suggests a general principle: any projection onto a basis can absorb the hardware's sampling grid into the filter function, a route worth testing for the other non-uniform grids used in geophysical models.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Quick take: this is a project deliverable that is more honest in the body than in the abstract. The GPU optimization work is real, detailed, and useful; the 23x speedup headline is not. The optical processor section is speculative but clearly labeled.

What's actually new: the report documents systematic GPU optimization of a spherical-harmonics spectral transform: batching variable-size cuBLAS DGEMMs with zero padding, masking to collapse an irregular loop nest with OpenACC collapse(3), flattening Fortran derived types so the compiler uses constant memory, and removing device synchronization between cuFFT calls. These are known techniques individually, but the combination and the practical pitfalls are described well enough to guide re-implementation. The MPDATA kernel result is cleaner: from 44GB/s to 344GB/s, 66% of STREAM, with a baseline that was clearly bandwidth-bound from the start. The CPU sections report modest, honest gains. The optical chapters are a genuine novelty -- a multichannel astigmatic optical processor -- and the authors are candid that quantitative extraction is work in progress.

Where it's soft: the headline 23x spectral transform speedup is measured against a GPU port that the report itself shows was dominated by Driver API calls and small kernels, with Compute sparse. The optimizations are exactly what you'd do to repair an immature first port: hoist allocations, expose all three parallel dimensions, batch library calls. So the 23x mostly measures the distance from a crippled port to a properly structured one; it does not indicate a 23x improvement in the spectral transform algorithm itself. The absolute 40%-of-roofline figure is more meaningful, but there are no raw timings, error bars, or repetitions, and the software artifact is only referenced as 'accessible via the ESCAPE Software Stash' with no URL. The abstract's 'demonstrated feasibility' for the optical processor overstates the chapter's own language: 4.4.7 calls high-fidelity numerical results 'challenging and still a work-in-progress,' and 4.4.13 withholds a key method for IP reasons. None of this is hidden, which I credit, but the summary should match the body.

Who it's for: someone working on GPU optimization of weather/climate codes, or evaluating alternative architectures for NWP. The MPDATA and GPU restructuring details earn a serious referee; the optical section is a proof-of-concept that needs substantially more validation. Recommend peer review with the condition that the authors supply artifacts and raw measurements, and re-baseline or reword the speedup claims.

Referee Report

4 major / 5 minor

Summary. This manuscript is the ESCAPE deliverable D3.3, reporting optimization results for weather and climate dwarfs on CPUs, Xeon Phi, GPUs, and the Optalysys optical processor. The main claims are: 1.1-2.5x speedups from compiler directives on CPU-based systems; 23x and 9x speedups on GPUs for the spherical-harmonics spectral transform and MPDATA, respectively; and demonstrated feasibility of a spectral transform in a purely optical fashion. The GPU work is placed in a roofline context, with the optimized spectral-transform time step reaching around 40% of the P100 roofline and the optimized MPDATA kernel reaching 66% of STREAM bandwidth. The optical section proposes two architectures, presents some experimental measurements, and discusses engineering hurdles and system-level caveats.

Significance. If the baseline-fairness concern is addressed, the report contains useful engineering data: the absolute roofline percentages, the description of loop restructuring, library batching, and constant-memory flattening give actionable guidance for GPU porting, and the CPU directive work is credible. The radiation refactoring results are also presented with an explicit acknowledgment that the huge speedups reflect legacy code, which is the right framing. The optical work is a genuine proof-of-concept with honest limitations, but the executive summary's 'demonstrated feasibility' language is stronger than the body supports. The main quantitative results are empirical measurements, not derived from a circular argument, and the later sections contain explicit caveats. The paper could be a valuable experience report if the headline claims are reframed and if raw timings and artifacts are made available.

major comments (4)

[§4.2.2, §4.2.3, §4.2.7, Executive Summary] The 23x speedup for the spectral transform is measured against a baseline that the report itself shows to be API-overhead-bound: Figure 2 is dominated by Driver API calls with sparse Compute activity, and Figure 1 places the original code 1-2 orders below the P100 roofline. The optimizations in §4.2.3, including hoisting allocations and copies out of the time-step loop, exposing all three levels of parallelism, batching cuBLAS/cuFFT calls, and removing per-call synchronization, are standard repairs for an immature first GPU port. The executive summary's formulation 'order of magnitude acceleration ... using data locality optimizations' therefore overstates what was demonstrated; the defensible quantitative result is the optimized code's roughly 40% roofline performance. The authors should either justify the original GPU port as a representative baseline or reframe the headline around absolute roofline performance and, where possible, compare against a reasonably optimized GPU implementation.
[§5.1.2, §5.1.4] The MPDATA speedup of 9.4x is also relative to a baseline that achieves only 44 GB/s, less than 10% of the STREAM bandwidth on the P100. The absolute result, 66% of STREAM bandwidth after optimization, is a good outcome, but it should be reported as a speedup over the initial OpenACC port rather than as an algorithmic improvement of the MPDATA scheme. As written, the Executive Summary's 9x figure inherits the same baseline-fairness concern as the spectral-transform number, even though the baseline here is less severely dominated by API overhead.
[§4.2.7, §5.1.4, §2.2.1] The headline numbers are not independently reproducible from the manuscript. No raw timings, repetitions, error bars, compiler and driver versions, or detailed platform settings are given for the GPU results; Figure 16 shows only relative bars. The text states that the software is accessible via the ESCAPE Software Stash, but no artifact identifier, repository URL, or access instructions are provided. Without these, the 23x, 9.4x, 40%-roofline, and 66%-bandwidth claims cannot be checked by a reader.
[§4.4.7, §4.4.2.2, §4.4.13, §4.4.14, Executive Summary] The executive summary's claim of 'demonstrated feasibility of a spectral transform in a purely optical fashion' is stronger than the evidence in the body. Section 4.4.7 states that obtaining high-fidelity numerical results is 'challenging and still a work-in-progress'; §4.4.2.2 explicitly reports a systematic error in the experimental measurements; §4.4.13 withholds the key sampling-compensation method for IP reasons; and §4.4.14.1 argues that a spectral-transform coprocessor 'does not necessarily offer a significant advantage' because the MPI_Alltoallv communication dominates in the IFS context. The optical section should be presented as a proof-of-concept with quantified accuracy limits, and the feasibility claim should be correspondingly qualified.

minor comments (5)

[§4.4.2.2, §4.4.8, Figure 30, Figure 21] There are several typos worth correcting: 'Legndre' in §4.4.2.2, 'implantation' in §4.4.8, 'iIntra' in the Figure 30 caption, and 'treads' in the Figure 21 caption.
[Figure 21 caption] The caption says 'The left graph represents all treads roofs and the left graph represents single thread roofs'; the second mention should almost certainly be 'right graph'.
[§4.3, §4.3.4] The subsection numbering repeats: '4.3.1 Summary' appears after §4.3.4, where it should be §4.3.5 or a similarly renumbered section.
[§6.3.4, Tables 3 and 4] Tables 3 and 4 are referenced in the text but their actual content appears to be missing from the submitted manuscript; either include the tables or remove the references.
[§10 References] The reference list contains only one entry, Juday 1993, while the text cites the ESCAPE website, the ECMWF dwarf documentation, an Intel tool, and Acosta's presentation by name or URL. A complete reference list should be added.

Circularity Check

0 steps flagged · score 0.0 of 10

No circular derivation: all claims are empirical measurements against external baselines and roofline ceilings.

full rationale

This deliverable reports measured speedups and roofline-relative performance; it contains no derivation chain in which an output quantity is defined in terms of the quantity it claims to predict. The GPU speedup of 23x for the spectral transform and 9.4x for the MPDATA kernel are direct before/after timings of specific code versions, not quantities fitted to the same data. The roofline comparisons use the externally defined STREAM benchmark as the bandwidth ceiling, and the paper explicitly states the measured 40% and 66% figures relative to that ceiling, so the ceiling is not constructed from the results. The optical section makes a feasibility claim supported by experimental demonstrations and explicitly hedged as work-in-progress in Sections 4.4.7 and 4.4.14, and its scaling arguments are conventional complexity comparisons rather than self-referential predictions. The references list only Juday (1993) for filter design, so no self-citation chain is load-bearing. The reader's concern that the baseline GPU port was immature and API-overhead-bound is a legitimate correctness or representativeness caveat about the comparison baseline, not circularity, because the baseline is an independently defined prior implementation rather than a fitted parameter or a renamed output. No step in the paper reduces, by its own equations or by citation, to its own inputs.

Assumptions & free parameters 2 free parameters · 3 assumptions · 1 invented entities

The central performance claims rely on standard engineering assumptions about measurement and hardware ceilings rather than on fitted constants. The only hand-tuned values are performance parameters (NPROMA, thread counts). The optical feasibility claim rests on the unverified assumption that measured intensity peaks correspond to accurate projections; the report itself flags systematic error.

free parameters (2)

NPROMA memory blocking size = 64 (on KNL 7250)
Chosen by parameter study in Section 8.1.3 to minimize execution time of the cloud microphysics dwarf; the reported speedup depends on this tuning.
OpenMP thread count = 128/134 on KNL 7250
Part of the same parameter study in Section 8.1.3; the reported cloud microphysics speedup depends on this execution configuration.

assumptions (3)

domain assumption The roofline model, with STREAM-measured bandwidth as the memory ceiling, provides a valid upper bound for achievable kernel performance.
Used throughout Sections 4.2 and 5.1 to argue that kernels are near their hardware limits.
domain assumption A lens performs an exact optical Fourier transform of the complex field at its back focal plane.
Foundational to the optical correlator design in Section 4.4.1.3.
domain assumption The central pixel of the optical correlation output equals the dot-product (projection) of input data with the basis function.
This is the load-bearing step for extracting Fourier-Legendre coefficients in Section 4.4.2.1; the report notes systematic error due to finite sampling.

invented entities (1)

Multichannel astigmatic optical processor
purpose: Cascades a 1D optical Fourier transform stage with a 2D optical correlator to evaluate parallel Legendre transforms for spherical harmonic decomposition.
Described as a prototype in Section 4.4.7 with qualitative output images; the report states 'quantitatively extracting different Fourier-Legendre coefficients is challenging' and no accuracy or throughput benchmark is provided.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Performance report and optimized implementation of Weather & Climate Dwarfs on GPU, MIC and Optalysys Optical Processor." pith.science (2026). https://pith.science/paper/KWY3PC6S

@misc{pith2026190806096,
  author       = {Pith},
  title        = {Pith review of: Performance report and optimized implementation of Weather & Climate Dwarfs on GPU, MIC and Optalysys Optical Processor},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/KWY3PC6S}},
  note         = {Machine review of arXiv:1908.06096}
}

read the original abstract

This document is one of the deliverable reports created for the ESCAPE project. ESCAPE stands for Energy-efficient Scalable Algorithms for Weather Prediction at Exascale. The project develops world-class, extreme-scale computing capabilities for European operational numerical weather prediction and future climate models. This is done by identifying Weather & Climate dwarfs which are key patterns in terms of computation and communication (in the spirit of the Berkeley dwarfs). These dwarfs are then optimised for different hardware architectures (single and multi-node) and alternative algorithms are explored. Performance portability is addressed through the use of domain specific languages. Here we summarize the work performed on optimizations of the dwarfs on CPUs, Xeon Phi, GPUs and on the Optalysys optical processor. We limit ourselves to a subset of the dwarf configurations and to problem sizes small enough to execute on a single node. Also, we use time-to-solution as the main performance metric. Multi-node optimizations of the dwarfs and energy-specific optimizations are beyond the scope of this report and will be described in Deliverable D3.4. To cover the important algorithmic motifs we picked dwarfs related to the dynamical core as well as column physics. Specifically, we focused on the formulation relevant to spectral codes like ECMWF's IFS code. The main findings of this report are: (a) Acceleration of 1.1x - 2.5x of the dwarfs on CPU based systems using compiler directives, (b) order of magnitude acceleration of the dwarfs on GPUs (23x for spectral transform, 9x for MPDATA) using data locality optimizations and (c) demonstrated feasibility of a spectral transform in a purely optical fashion.

Figures

Figures reproduced from arXiv: 1908.06096 by the authors.

**Figure 5.** The performance of the SH dwarf TL159 test case on the NVIDIA Tesla P100 GPU, in the context of [PITH_FULL_IMAGE:figures/full_fig_p005_5.png] view at source ↗

**Figure 16.** The time taken by the original and optimized GPU versions of the compute_fluxzdiv kernel on the [PITH_FULL_IMAGE:figures/full_fig_p006_16.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

15 extracted references · 15 canonical work pages

[1]

Through characterisation of the devices, and compensating for deviations during operation of the device

work page
[2]

Profiling and Computational Performance of IFS using BSC Tools

Through in situ optimisation of the device. Corrections are represented agnostically by some appropriate parameterisation, and these parameters are established through some in situ process. Optalysys has experience with the former technique. While it does represent a ‘gold standard’ – in that the error that is being compensated for can be explicitly measu...

work page 2017
[3]

do jnode = 1,this%geom%nb_nodes

work page
[4]

do jlev = 1,this%geom%nb_levels

work page
[5]

do jedge = 1,inode2edges_size(jnode)

work page
[6]

iedge = inode2edges(jnode, jedge)

work page
[7]

zadd = real(this%geom%node2edge_sign(jedge,jnode),wp)

work page
[8]

zsum = zsum+zadd*pFx(jlev,iedge)

work page

Show all 15 references

[9]

pdivVD(jlev,jnode) = zsum/pvol(jnode) &
[10]

& +(pFz(jlev+1,jnode)-pFz(jlev,jnode))/this%dz
[12]

roofline

enddo It can clearly be seen that there is a relatively low arithmetic intensity (where for more details on this please see the Spherical Harmonics GPU section), so global memory bandwidth will be the limiting factor for test cases of realistic size and the performance “roofli...

2017
[13]

X: portable data-structures, tuned for Intel Xeon/Xeon Phi
[14]

G: portable data-structures, tuned for GPU
[15]

GNM: non-portable data-structures, tuned solely for the GPU This is all described in more details in the paper. 6.3.2 Performance Results GPU P100 timings of the G and GNM codes on the 400x400x80 test case were 4s and 1.7s, respectively, which perfectly matches our naive expec...

2017
[2016]

nabla operators divergence 2d

This is shown in Table 3. ESCAPE 2017 D3.3 Performance report and optimized implementation of Weather & Climate Dwarfs on GPU, MIC and Optalysys Optical Processor 55 Table 4 Performance improvements by individual node refactoring. Figure 17 Impact of different data layouts on ...

1993

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.

[1] [1]

Through characterisation of the devices, and compensating for deviations during operation of the device

work page

[2] [2]

Profiling and Computational Performance of IFS using BSC Tools

Through in situ optimisation of the device. Corrections are represented agnostically by some appropriate parameterisation, and these parameters are established through some in situ process. Optalysys has experience with the former technique. While it does represent a ‘gold standard’ – in that the error that is being compensated for can be explicitly measu...

work page 2017

[3] [3]

do jnode = 1,this%geom%nb_nodes

work page

[4] [4]

do jlev = 1,this%geom%nb_levels

work page

[5] [5]

do jedge = 1,inode2edges_size(jnode)

work page

[6] [6]

iedge = inode2edges(jnode, jedge)

work page

[7] [7]

zadd = real(this%geom%node2edge_sign(jedge,jnode),wp)

work page

[8] [8]

zsum = zsum+zadd*pFx(jlev,iedge)

work page