Scientific Computing in the AI Era

The 2019 edition of this manual taught shebang lines, OpenMP pragmas on CPU loops, MPI in C and Fortran, Conda and Jupyter. All of that still works, and most of it still matters. What changed is where the flops live. The large systems are built around GPUs and other accelerators, the codes that matter most have been rewritten or re-plumbed to run on them, and a new class of scientific tools (surrogate models, neural weather forecasters, protein structure predictors) now sits next to the simulation codes and competes with them for the same nodes. This chapter covers the working scientist's and facilitator's view: what classic simulation looks like now, which GPU programming model to pick, how MPI has evolved, how to port a code without wasting a year, which applications are ready, and how AI for science fits in. Hardware details live in Accelerators and Interconnects; the software stack underneath (drivers, CUDA/ROCm, containers) is in Systems Software.

Where Classic Simulation Stands

Computational fluid dynamics, climate and weather, molecular dynamics, materials and electronic structure, and astrophysics still consume a large share of research computing cycles. The equations have not changed. What changed is the hardware those equations are solved on, and the pressure to use AI models as a complement to (or replacement for) some runs.

The June 2026 TOP500 list is a useful snapshot. The new No. 1, LineShine in Shenzhen, is a CPU-only design at 2.198 Exaflop/s HPL, a reminder that CPU machines are not dead. The next four, El Capitan (AMD Instinct MI300A accelerators, 1.809 Exaflop/s), Frontier (1.353 Exaflop/s), Aurora (Intel Data Center GPU Max accelerators, 1.012 Exaflop/s) and JUPITER Booster, are accelerator machines.[1] The difference shows in mixed precision: LineShine's HPL-MxP result is 3.6x its HPL score, while the accelerated systems run 6.3x to 11.5x.[1][2] For most research centers, though, the practical reality is that new capacity arrives as GPU nodes, and codes that cannot use them are left on shrinking CPU partitions.

What maps well to GPUs

• Bandwidth-bound stencils and sparse linear algebra (CFD, climate dynamical cores, lattice methods). These are limited by memory bandwidth, and HBM gives GPUs several times the bandwidth of a CPU socket. See Memory and Storage.
• Dense linear algebra and FFTs (electronic structure, quantum chemistry, PME electrostatics in MD). These map onto vendor libraries (cuBLAS, cuFFT, rocBLAS, rocFFT, oneMKL) that are already tuned.
• Particle methods with lots of independent work (MD non-bonded forces, N-body, Monte Carlo transport).

What still struggles

• Codes with deep, branchy, pointer-chasing logic per element (some chemistry and physics parameterizations, adaptive mesh bookkeeping).
• Small problems that cannot fill a modern GPU. A 20,000-atom system or a coarse 2D mesh may run faster on a well-used CPU node than on a mostly idle GPU.
• Workflows dominated by I/O or by serial pre- and post-processing, where the GPU sits idle most of the job.
• Codes that need FP64 everywhere on hardware whose FP64 rate has not kept up with its low-precision tensor throughput. See Mixed Precision below and The Convergence.

GPU Programming Models

There is no single right answer, and the choice outlives any one machine. The honest decision criteria are: which vendors you must support in five years, what language the code is already in, how much of the team's time goes to performance engineering, and whether the code has a large community that has already chosen for you.

Native models: CUDA, HIP and SYCL

CUDA is still the reference point. Most GPU libraries, profilers, tutorials and Stack Overflow answers assume it. If a code only ever needs to run on NVIDIA hardware, CUDA plus the CUDA-X libraries is the lowest-risk path to top performance.

HIP is AMD's C++ runtime API and kernel language, deliberately close to CUDA. AMD's HIPIFY tools (hipify-clang, a Clang-based translator, and hipify-perl, a text-based one) convert CUDA sources to HIP as a starting point.[3] Many projects keep one source tree and a thin macro layer that maps to either CUDA or HIP.

SYCL is a Khronos standard for heterogeneous programming in ISO C++. Intel's oneAPI DPC++/C++ compiler fully supports SYCL 2020, and AdaptiveCpp is an independent implementation that targets multiple vendors.[4] GROMACS, for example, recommends SYCL builds for both AMD (with AdaptiveCpp) and Intel GPUs.[5]

Directive models: OpenMP offload and OpenACC

Directives let you keep the loop you already have and annotate it. That is attractive for large Fortran codes with decades of validated physics. OpenMP 6.0 was released at SC24 in November 2024; its device-side additions include a workdistribute directive for Fortran array syntax, structured asynchronous data-mapping regions, better control of allocatable variables and per-team device memory via groupprivate.[6] OpenACC 3.4 was announced at ISC 2025; its most visible addition is a capture modifier on data clauses so that arrays can still be placed in discrete device memory on shared or unified memory systems, fixing race conditions seen in real codes.[7]

In practice, OpenACC is best supported by the NVIDIA HPC SDK compilers, while OpenMP offload is the vendor-neutral directive path (NVIDIA, AMD, Intel, HPE Cray and LLVM compilers all implement it, with varying completeness). VASP is a good illustration: it ships an OpenACC port for NVIDIA GPUs and an OpenMP offload port for AMD and Intel GPUs.[8]

Portability layers: Kokkos and RAJA

Kokkos is a C++ library that lets one kernel implementation run on different hardware through compile-time backends, and it can also adjust data layout per architecture. Kokkos 5.0 (November 2025) moved the minimum standard to C++20, and the project ships on a roughly four-month cadence (5.2 arrived in July 2026).[9][10] LAMMPS and the E3SM atmosphere model are two large users.[11][12] RAJA, from LLNL, is a comparable set of C++ loop abstractions for architecture and programming-model portability.[13] Choose one of these when a C++ code must run well on NVIDIA, AMD and Intel GPUs and on CPUs, and the team is prepared to adopt its idioms.

Standard language parallelism

C++17 parallel algorithms (std::transform, std::for_each with an execution policy) and Fortran do concurrent are part of the languages themselves. With the NVIDIA HPC SDK, the -stdpar option (default -stdpar=gpu) offloads these constructs to NVIDIA GPUs; GPU acceleration of do concurrent arrived in nvfortran with HPC SDK 20.11.[14][15] The appeal is obvious: no directives, no new language, and the code still compiles and runs on a CPU with any conforming compiler. The limits are real too: you get less control over data placement and kernel shape, and support on non-NVIDIA compilers is less mature.

Julia and Python

Julia has first-class GPU packages: CUDA.jl for NVIDIA, AMDGPU.jl for ROCm (described by the JuliaGPU project as somewhat less developed), oneAPI.jl for Intel GPUs and Metal.jl for Apple Silicon, the last two with more limited library integration.[16]

In Python, the options have sorted themselves out:

• CuPy is a NumPy/SciPy-compatible array library for NVIDIA CUDA and AMD ROCm; on ROCm it maps CUDA library calls to their HIP equivalents.[17] It is the fastest way to move array-heavy analysis to a GPU.
• Numba compiles Python functions, but its built-in CUDA target is deprecated, with development moved to NVIDIA's separate numba-cuda package.[18] New GPU kernels written with Numba should target numba-cuda.
• JAX traces Python functions and hands them to the XLA compiler, which fuses operations and specializes for CPU, GPU or TPU.[19] Automatic differentiation plus XLA make it popular for differentiable simulation and for AI-for-science models such as GraphCast and GenCast.
• cuPyNumeric (NumPy on the Legate runtime, scaling from one GPU to many nodes) has reached end of life. The final release is v26.06.01 and no further releases, fixes or support are planned.[20] Do not start new projects on it; existing users should plan a migration to CuPy, JAX or Dask-based alternatives.

Field note: The programming model decision is really a staffing decision. A directive port can be maintained by domain scientists who mostly write physics. A CUDA or Kokkos code needs at least one person who is comfortable reading profiler output and rewriting kernels. Pick the model your group can keep alive after the postdoc who did the port leaves.

MPI Today

MPI remains the way nearly every multi-node simulation code communicates, and it has kept evolving. The MPI Forum approved MPI 4.1 on November 2, 2023 and MPI 5.0 on June 5, 2025; work on a follow-on ("MPI Next") with features still undefined is ongoing.[21][22] MPI 4.1 was scoped as minor additions and corrections plus readability work.[23] MPI 5.0 is the bigger change: it adds a standard Application Binary Interface (ABI), defining the layout of the status object, the type of handles such as MPI_Comm, literal values for many constants, and new functions MPI_ABI_GET_VERSION and MPI_ABI_GET_INFO, plus a Fortran ABI section.[24]

Why the ABI matters to a research computing center: today an application built against Open MPI will not run against MPICH without a rebuild, which multiplies the software stack and complicates containers. A standard ABI means a binary (or a container) built once can, in principle, run against whichever conforming MPI library the site provides. MPICH 5.0.0 is the first stable MPICH release series with full MPI-5 support including the ABI.[25] Expect other implementations to follow, and expect a transition period where you still need implementation-specific builds.

GPU-aware MPI

A GPU-aware (CUDA-aware, ROCm-aware) MPI accepts device pointers directly in send, receive and collective calls. Without it, every halo exchange means copying data from GPU to host, calling MPI, and copying back. With it, the library can use GPUDirect RDMA between the NIC and GPU memory across nodes and peer-to-peer IPC within a node.

• Open MPI 5.0.x supports CUDA buffers either through UCX (the preferred mechanism, which requires UCX itself built with CUDA) or through its internal CUDA support; it is configured with --with-cuda, detects CUDA components at run time, and exposes compile-time and run-time checks (MPIX_CUDA_AWARE_SUPPORT, MPIX_Query_cuda_support() in mpi-ext.h).[26] ROCm device buffers are supported through UCX, which is the only formally supported transport for ROCm in the 5.0.x series.[27]
• MPICH enables GPU support automatically when CUDA, Level Zero (ZE) or HIP runtimes are found at configure time (--with-cuda, --with-ze, --with-hip), and it can be disabled at run time with MPIR_CVAR_ENABLE_GPU=0 to avoid GPU initialization overhead in CPU-only jobs.[28]
• HPE Cray MPICH supports NVIDIA and AMD GPUs using GPU-NIC RDMA between nodes and GPU peer-to-peer IPC within a node, but you must set MPICH_GPU_SUPPORT_ENABLED=1 at run time.[29] Forgetting it is the single most common cause of crashes when a code passes device pointers on a Cray EX system.

A minimal launch on a Slurm GPU cluster looks like the following. Partition names, module names and counts are examples.

#!/bin/bash
#SBATCH --job-name=halo-test
#SBATCH --partition=gpu            # example partition name
#SBATCH --nodes=2
#SBATCH --ntasks-per-node=4        # one MPI rank per GPU
#SBATCH --gpus-per-node=4
#SBATCH --cpus-per-task=8
#SBATCH --time=00:30:00

module load cuda openmpi           # example module names; use your site's GPU-aware MPI

# On HPE Cray systems with cray-mpich, this is required for device buffers:
# export MPICH_GPU_SUPPORT_ENABLED=1

srun --gpu-bind=closest ./halo_exchange

Scheduling details (GPU binding, MPS, MIG) are in Scheduling and Orchestration. How collectives behave on the fabric, and how NCCL and RCCL relate to MPI, is in Interconnects.

Porting a Code to GPUs

The failures we see in porting projects are rarely about syntax. They come from porting the wrong thing, or from moving data back and forth so often that the GPU never gets a chance to help. A workable sequence:

1. Build a regression test first. Pick two or three representative inputs and define what "same answer" means (bitwise is usually impossible; agree on tolerances with the scientists).
2. Profile the CPU code. Find the handful of routines that take most of the time. If 80 percent of runtime is in three loops, port those; do not touch the input parser.
3. Make data resident. Allocate the main arrays on the GPU once, keep them there across time steps, and only copy what the host truly needs (diagnostics, I/O, halos if MPI is not GPU-aware).
4. Port kernels, then fuse. Get correctness with simple one-to-one kernels, then merge small kernels that touch the same data.
5. Move communication onto the device. Pack halos on the GPU and pass device pointers to a GPU-aware MPI.
6. Profile again on the GPU. NVIDIA Nsight Systems (for example nsys profile --stats=true ./app) shows a timeline of kernels, copies and MPI calls; gaps and long host-device copies jump out immediately.[30] AMD's ROCm profilers and Intel's VTune fill the same role on their hardware. More on this in Observability.

Field note: The first GPU build of a legacy code is often slower than the CPU version. That is normal and almost always caused by host-device transfers inside the time loop. Look at the profiler timeline before touching a single kernel. If you see copies between every kernel, fix data residency first.

SAXPY in CUDA and in OpenMP offload

SAXPY (y = a*x + y) is the "hello world" of GPU computing. In CUDA, you write the kernel explicitly and manage device memory. This version uses unified (managed) memory to keep the example short; production codes usually manage device allocations explicitly.

// saxpy.cu -- build: nvcc -O3 -arch=sm_90 saxpy.cu -o saxpy   (sm_90 = Hopper; example)
#include <cstdio>
#include <cuda_runtime.h>

__global__ void saxpy(int n, float a, const float *x, float *y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) y[i] = a * x[i] + y[i];
}

int main() {
    const int n = 1 << 24;
    float *x, *y;
    cudaMallocManaged(&x, n * sizeof(float));
    cudaMallocManaged(&y, n * sizeof(float));
    for (int i = 0; i < n; ++i) { x[i] = 1.0f; y[i] = 2.0f; }

    int threads = 256;
    int blocks = (n + threads - 1) / threads;
    saxpy<<<blocks, threads>>>(n, 3.0f, x, y);
    cudaDeviceSynchronize();

    printf("y[0] = %f (expect 5.0)\n", y[0]);
    cudaFree(x); cudaFree(y);
    return 0;
}

The same computation with OpenMP offload is a single directive on the existing loop. The map clauses say what moves to and from the device.

/* saxpy_omp.c
 * NVIDIA HPC SDK:  nvc -O3 -mp=gpu saxpy_omp.c -o saxpy_omp
 * LLVM clang:      clang -O3 -fopenmp --offload-arch=sm_90 saxpy_omp.c -o saxpy_omp  (example arch)
 */
#include <stdio.h>
#include <stdlib.h>

int main(void) {
    const int n = 1 << 24;
    float a = 3.0f;
    float *x = malloc(n * sizeof(float));
    float *y = malloc(n * sizeof(float));
    for (int i = 0; i < n; ++i) { x[i] = 1.0f; y[i] = 2.0f; }

    #pragma omp target teams distribute parallel for map(to: x[0:n]) map(tofrom: y[0:n])
    for (int i = 0; i < n; ++i)
        y[i] = a * x[i] + y[i];

    printf("y[0] = %f (expect 5.0)\n", y[0]);
    free(x); free(y);
    return 0;
}

The -mp=gpu flag is how the NVIDIA compilers offload OpenMP regions, and LLVM uses --offload-arch to name the target GPU (for example sm_80 for NVIDIA A100 or gfx90a for AMD MI200-series).[14][31] In Fortran with standard parallelism, the loop needs no directive at all:

! saxpy_dc.f90 -- build: nvfortran -O3 -stdpar=gpu saxpy_dc.f90 -o saxpy_dc
program saxpy_dc
  implicit none
  integer, parameter :: n = 2**24
  real, allocatable :: x(:), y(:)
  real :: a
  integer :: i
  allocate(x(n), y(n))
  a = 3.0
  x = 1.0
  y = 2.0
  do concurrent (i = 1:n)
     y(i) = a * x(i) + y(i)
  end do
  print *, 'y(1) =', y(1), ' (expect 5.0)'
end program saxpy_dc

A stencil with resident data

Real codes iterate. The key idea is a target data region that keeps arrays on the device for the whole solve, so the loop body does no host-device copies. This is a 2D Jacobi relaxation, the core of many CFD and diffusion solvers.

/* One Jacobi sweep per iteration; u and unew stay on the GPU the whole time. */
void jacobi(int nx, int ny, int iters, double *restrict u, double *restrict unew) {
    #pragma omp target data map(tofrom: u[0:nx*ny]) map(alloc: unew[0:nx*ny])
    {
        for (int it = 0; it < iters; ++it) {
            #pragma omp target teams distribute parallel for collapse(2)
            for (int j = 1; j < ny - 1; ++j)
                for (int i = 1; i < nx - 1; ++i)
                    unew[j*nx + i] = 0.25 * (u[j*nx + i - 1] + u[j*nx + i + 1]
                                           + u[(j-1)*nx + i] + u[(j+1)*nx + i]);

            #pragma omp target teams distribute parallel for collapse(2)
            for (int j = 1; j < ny - 1; ++j)
                for (int i = 1; i < nx - 1; ++i)
                    u[j*nx + i] = unew[j*nx + i];
        }
    }   /* u is copied back to the host once, here */
}

A GPU-aware MPI halo exchange

Once the stencil is split across ranks (one per GPU), each rank must exchange boundary rows with its neighbors every iteration. With a GPU-aware MPI, you hand MPI the device addresses. In OpenMP, a use_device_ptr clause on a target data construct replaces a host pointer with the corresponding device pointer inside the region, so host code (here, MPI calls) sees device addresses. This sketch assumes a 1D decomposition in j with one ghost row on each side, and that u is already mapped by an enclosing target data region.

#include <mpi.h>

/* Exchange ghost rows with up/down neighbors. u has (ny_local + 2) rows of nx doubles.
 * up/down are neighbor ranks or MPI_PROC_NULL at the domain edges. */
void halo_exchange(double *u, int nx, int ny_local, int up, int down, MPI_Comm comm) {
    MPI_Request req[4];
    #pragma omp target data use_device_ptr(u)
    {
        /* Inside this region, u refers to the DEVICE copy, so MPI gets device pointers. */
        double *first_real = u + 1 * nx;          /* row 1          */
        double *last_real  = u + ny_local * nx;   /* row ny_local   */
        double *ghost_lo   = u;                   /* row 0          */
        double *ghost_hi   = u + (ny_local + 1) * nx;

        MPI_Irecv(ghost_lo,   nx, MPI_DOUBLE, down, 0, comm, &req[0]);
        MPI_Irecv(ghost_hi,   nx, MPI_DOUBLE, up,   1, comm, &req[1]);
        MPI_Isend(first_real, nx, MPI_DOUBLE, down, 1, comm, &req[2]);
        MPI_Isend(last_real,  nx, MPI_DOUBLE, up,   0, comm, &req[3]);
        MPI_Waitall(4, req, MPI_STATUSES_IGNORE);
    }
}

Contiguous rows are the easy case. For column halos in a 2D decomposition, pack the strided data into a contiguous device buffer with a small kernel first, send that, and unpack on the receiving side. Overlap interior computation with the exchange once correctness is established. If the MPI library is not GPU-aware, this code will fail or silently read garbage; check with the implementation's query mechanism before trusting it.

GPU-Accelerated Applications

Most users never port anything. They run community codes, and the question is whether their code and their particular feature set run well on GPUs. The status of major applications, from the projects' own documentation:

Field note: "Runs on GPUs" is feature-specific. A code may accelerate the standard path beautifully and fall back to the CPU for the one thermostat, force field or solver option your group depends on. Before you promise a lab that GPUs will help, run their actual input deck on one GPU node and one CPU node and compare time to solution and cost per result.

AI for Science

AI for science is not one thing. It is at least four different workloads, each with different compute and validation needs.

Surrogates, physics-informed ML and neural operators

A surrogate is a model trained on simulation (or experimental) output that predicts results far faster than rerunning the simulation. It is useful for design sweeps, uncertainty quantification and inverse problems, and it is only trustworthy inside the range of its training data.

Physics-informed neural networks (PINNs), introduced by Raissi, Perdikaris and Karniadakis in the Journal of Computational Physics in 2019, embed the governing PDE into the loss function so the network is penalized for violating the physics, and they handle both forward and inverse problems.[42] They are appealing when data are scarce but can be hard to train for stiff or multiscale problems.

Neural operators learn a map between function spaces (for example from initial conditions to the solution field), rather than a single solution. The Fourier Neural Operator paper reported modeling turbulent flow with zero-shot super-resolution and speeds up to three orders of magnitude faster than traditional PDE solvers on its test problems.[43] NVIDIA PhysicsNeMo packages these families (neural operators, GNNs, transformers, PINNs and hybrids) as PyTorch modules with scalable training pipelines.[41]

The compute pattern matters for capacity planning: generating training data is a classic HPC simulation campaign (often FP64, MPI, lots of storage), training is a GPU deep learning job (see Training at Scale), and inference is cheap enough to run interactively. One project can need all three in the same month.

AI weather and climate models

Weather is where machine-learned models moved fastest from papers into operations.

Two practical lessons. First, these models are trained on reanalysis data produced by physics models, so they depend on the physics pipeline rather than replacing it. ECMWF explicitly runs AIFS side by side with IFS. Second, inference is cheap but training is not: training runs use many GPUs or TPUs for weeks, while a lab can run inference on one GPU node. Most university groups will fine-tune or run these models, not train them from scratch.

Protein structure and biomolecular models

AlphaFold 2 (Nature, 2021) changed structural biology, and AlphaFold 3 (Nature, 2024) extended prediction to complexes of proteins with nucleic acids, small molecules and ions.[49][50] For research computing staff, licensing and workflow are the operational issues. The AlphaFold 3 inference code is on GitHub under Apache 2.0, but the model parameters are only available directly from Google under separate terms of use, so each user or lab must obtain them and a center cannot simply redistribute a ready-to-run container with weights.[51]

Open alternatives have matured. OpenFold3-preview, from the AlQuraishi Lab at Columbia and the OpenFold consortium, aims to be a bitwise reproduction of AlphaFold 3 and is Apache 2.0 licensed, with training data released.[52] Boltz-2 (MIT license, model weights and training pipeline open) jointly predicts structure and binding affinity; its developers report affinity accuracy approaching physics-based free-energy perturbation (FEP) at more than 1000x the speed.[53] Treat that as a claim to validate on your own targets, not a reason to retire FEP.

These tools have an awkward shape for schedulers: the multiple sequence alignment step is CPU and I/O heavy against multi-terabyte databases, and the model inference step wants one GPU with a lot of memory. Split them into separate job steps or a workflow with dependencies so GPUs are not held idle during database searches. Stage the databases on fast local or parallel storage.

Scientific foundation models and LLM agents

Aurora is one example of a broader trend toward large pretrained models for a scientific domain that are then fine-tuned for specific tasks. Similar efforts exist for materials, chemistry, genomics and climate. At the other end, general-purpose LLMs are being wrapped into research agents. Google's Co-Scientist, for example, is a multi-agent system built on Gemini that generates, critiques and refines hypotheses using extra test-time compute.[54] In the US, the Genesis Mission, launched by executive order on November 24, 2025 and led by DOE, aims to apply the national laboratories' supercomputers and data resources to AI-driven science.[55]

For a campus center, the immediate implications are concrete: researchers will ask for GPU time to fine-tune domain models, for inference endpoints they can call from notebooks and workflows (see Inference at Scale), and for agent frameworks that submit jobs on their behalf. The last one needs policy before it needs hardware: an agent with a user's credentials can burn an allocation or touch data that user can reach. Cover this in your acceptable use policy (see Operations).

Mixed Precision and FP64 Emulation

Modern accelerators devote most of their silicon to low-precision tensor math. Scientific codes that assume FP64 everywhere leave much of that capacity unused. Two approaches are moving into the mainstream.

Mixed-precision algorithms do most of the work in lower precision and recover full accuracy with a small amount of high-precision work, for example iterative refinement in linear solvers. The HPL-MxP benchmark measures exactly this. On its results list, El Capitan reports 16.68 Exaflop/s mixed precision against 1.809 Exaflop/s on FP64 HPL, a 9.2x ratio.[2] Real applications see far smaller gains than a benchmark designed for the purpose, but the direction is clear.

Floating-point emulation computes FP32 or FP64 matrix products by splitting each value into several lower-precision pieces, multiplying those on tensor cores, and recombining. NVIDIA introduced emulation in cuBLAS with CUDA 12.9 (FP32 emulated with BF16, called BF16x9), and CUDA 13.0 Update 2 added FP64 emulation on Blackwell INT8 tensor cores using Ozaki-scheme methods, with an automatic mode that picks parameters to match or beat native FP64 accuracy and falls back to native FP64 when emulation offers no benefit.[56][57] NVIDIA reports a 1.5x end-to-end speedup for Quantum ESPRESSO with the automatic setting on an RTX PRO 6000 Blackwell Server Edition, with results consistent to 12 significant digits, and a 2.4x matrix-product speedup in the ecTrans weather library with FP32 emulation.[57]

What this means in practice: emulation only helps GEMM-heavy codes that call cuBLAS (directly or through a library), it helps most on GPUs with weak native FP64, and it changes numerical behavior enough that you must rerun your validation suite. It does nothing for stencils, sparse solvers or particle codes that are bandwidth-bound. When evaluating hardware for FP64 science, ask vendors for native FP64 figures and for emulated results on your own workloads, and treat the two separately. Precision formats are covered in The Convergence.

Reproducibility

GPU codes make reproducibility harder in two ways. Numerically, parallel reductions and atomic operations change the order of floating-point additions from run to run, and GPU math can differ from CPU math. The NAMD release notes, for example, warn that GPU forces differ slightly from CPU-only calculations, most visibly in reported pressure.[34] Environmentally, a result now depends on the driver, the CUDA or ROCm version, the MPI build, and library versions, not just the source code.

• Define correctness as tolerance, not bits. Agree with the science team on which outputs must match and how closely, and record it with the regression tests.
• Pin the software stack. Spack environments capture which packages to build, how they are configured and the concretized result in spack.yaml and spack.lock, so a stack can be rebuilt exactly.[58]
• Use containers for portability, carefully. Apptainer's --nv and --rocm flags bring the host's GPU driver libraries into the container.[59] The container holds the user-space stack, but the host driver must still be new enough for the CUDA or ROCm version inside, and MPI inside the container must be compatible with the host fabric (the MPI 5.0 ABI should eventually ease this).
• Record provenance automatically. Log module lists, container digests, GPU model, driver version and job IDs with every production run. Workflow tools and simple wrapper scripts both work; what matters is that nobody has to remember to do it.
• For AI models, version the weights and the data. A surrogate or structure predictor is part of the method. Record the exact checkpoint, its license, and the training data version.

Interactive Computing at Scale: Jupyter, Dask and Ray

Notebooks became the front door for a large share of users. The questions for a center are how to give people notebooks on real compute without letting them camp on GPUs, and how to let Python scale past one node.

Jupyter on clusters

Open OnDemand provides JupyterLab as an interactive app: the user fills in a form, the portal submits a batch job, and the notebook runs on a compute node inside a normal Slurm allocation.[60] That keeps notebooks under the scheduler's accounting and limits. Pair it with short default wall times on GPU partitions, idle culling, and a small shared or fractional GPU partition (MIG or MPS; see Scheduling) for development, so full GPUs go to batch work.

Dask and Ray

Dask scales NumPy- and pandas-style work across many workers. The dask-jobqueue project deploys Dask on PBS, Slurm, MOAB, SGE, LSF and HTCondor, either as dynamic clusters where workers are submitted as separate jobs or as batch runners inside one allocation.[61]

# Dask on Slurm from a notebook or login node (values are examples)
from dask_jobqueue import SLURMCluster
from dask.distributed import Client

cluster = SLURMCluster(
    queue="standard",            # Slurm partition
    account="mylab",             # Slurm account
    cores=16,
    memory="64GB",
    walltime="02:00:00",
    job_extra_directives=["--constraint=intel"],   # example extra sbatch line
)
cluster.scale(jobs=4)            # submit 4 worker jobs
client = Client(cluster)

import dask.array as da
x = da.random.random((200_000, 20_000), chunks=(10_000, 20_000))
print(x.mean().compute())

Ray targets distributed Python more generally (task and actor parallelism, hyperparameter tuning, serving). Ray expects a head node plus workers, which does not match Slurm's model naturally; Ray 2.49 and later add a ray symmetric-run command that starts the Ray cluster on all allocated nodes and runs the entry script only on the head node, which simplifies running Ray inside an sbatch job.[62]

Field note: Dynamic Dask clusters submit many small jobs, which can look like abuse to a busy scheduler and can stall when the queue is long. For large runs, prefer the batch-runner pattern (one allocation, Dask inside it). Reserve dynamic clusters for interactive exploration on partitions with short waits.

Key Takeaways

• New capacity arrives mostly as GPU nodes. Bandwidth-bound and dense-math codes benefit most; small or branchy workloads may still belong on CPUs.
• Choose a GPU programming model by portability needs and staffing: CUDA for NVIDIA-only performance, OpenMP offload or OpenACC for incremental Fortran ports, Kokkos or RAJA for portable C++, standard parallelism for simple loops. cuPyNumeric is end of life; use CuPy, JAX or Dask instead.
• MPI 5.0 (June 2025) adds a standard ABI, which should eventually let one binary or container run against any conforming MPI. MPICH 5.0 already implements it.
• GPU-aware MPI is essential for multi-GPU codes. Verify it is built and enabled (for Cray MPICH, set MPICH_GPU_SUPPORT_ENABLED=1).
• Porting is mostly about data residency: profile first, keep arrays on the device, and move halo exchanges onto the GPU.
• AI weather and structure-prediction models are now operational tools, but they are trained on physics-derived data and need validation; check model weight licenses before building shared services.
• FP64 emulation in cuBLAS can speed up GEMM-heavy codes on hardware with weak native FP64, but it requires revalidation and does nothing for bandwidth-bound kernels.
• Reproducibility now means pinning drivers, libraries, containers and model checkpoints, and defining correctness by tolerance.

Sources

1. TOP500 June 2026 list and highlights (TOP500, 2026)
2. HPL-MxP results (HPL-MxP, 2026)
3. HIPIFY documentation (AMD ROCm, 2026)
4. SYCL overview and implementations (Khronos Group, 2026)
5. GROMACS 2025.0 installation guide: GPU acceleration options (GROMACS, 2025)
6. OpenMP ARB releases OpenMP 6.0 (OpenMP ARB, 2024)
7. Announcing OpenACC 3.4 at ISC 2025 (OpenACC, 2025)
8. GPU ports of VASP (VASP Wiki, 2026)
9. Kokkos 5.0: Embracing C++20 (High Performance Software Foundation, 2025)
10. Kokkos releases (Kokkos project, 2026)
11. KOKKOS package (LAMMPS documentation, 2026)
12. Modeling at Exascale: EAMxx (E3SM project, DOE, 2025)
13. RAJA documentation (Lawrence Livermore National Laboratory, 2026)
14. NVIDIA HPC Compilers User's Guide (NVIDIA, 2026)
15. Accelerating Fortran DO CONCURRENT with GPUs and the NVIDIA HPC SDK (NVIDIA, 2020)
16. JuliaGPU: GPU programming in Julia (JuliaGPU, 2026)
17. CuPy overview (CuPy documentation, 2026)
18. Numba for CUDA GPUs: built-in CUDA target deprecation notice (Numba documentation, 2026)
19. Just-in-time compilation (JAX documentation, 2026)
20. cuPyNumeric repository: end-of-life notice (NVIDIA, 2026)
21. MPI Documents (MPI Forum, 2025)
22. MPI Forum home: current efforts (MPI Forum, 2025)
23. MPI 4.1 scope (MPI Forum, 2023)
24. MPI 5.0 standard: Changes in MPI-5.0 (MPI Forum, 2025)
25. MPICH v5.0.0 release notes (MPICH project, 2026)
26. Open MPI 5.0.x: CUDA (Open MPI documentation, 2026)
27. Open MPI 5.0.x: ROCm (Open MPI documentation, 2026)
28. MPICH 4.3.0 README: GPU support (MPICH project, 2025)
29. intro_mpi: GPU support in Cray MPICH (HPE Cray Programming Environment 25.09, 2025)
30. Nsight Systems User Guide (NVIDIA, 2026)
31. OpenMP command-line argument reference (LLVM OpenMP, 2026)
32. Heterogeneous parallelization and GPU acceleration (GROMACS, 2026)
33. Getting good performance from mdrun (GROMACS documentation, 2026)
34. NAMD 3.0 release notes (University of Illinois, 2024)
35. The pmemd.cuda GPU implementation (Amber, 2026)
36. Quantum ESPRESSO: One Further Step toward the Exascale (J. Chem. Theory Comput., 2023)
37. Building an Accelerated OpenFOAM Proof-of-Concept Application using Modern C++ (arXiv, 2025)
38. GPU Acceleration project (NCAR, 2026)
39. AceCAST documentation: GPU-accelerated WRF (TempoQuest, 2026)
40. Optimizing the Weather Research and Forecasting Model with OpenMP Offload and Codee (arXiv, 2024)
41. NVIDIA PhysicsNeMo overview (NVIDIA, 2026)
42. Raissi, Perdikaris, Karniadakis: Physics-informed neural networks (Journal of Computational Physics 378, 2019)
43. Li et al.: Fourier Neural Operator for Parametric Partial Differential Equations (arXiv, 2020)
44. Lam et al.: GraphCast: Learning skillful medium-range global weather forecasting (arXiv / Science, 2023)
45. Price et al.: Probabilistic weather forecasting with machine learning (Nature, 2025)
46. Bodnar et al.: A Foundation Model for the Earth System (Nature 641, 2025)
47. FourCastNet 3: A geometric approach to probabilistic machine-learning weather forecasting at scale (arXiv, 2025)
48. ECMWF's ensemble AI forecasts become operational (ECMWF, 2025)
49. Jumper et al.: Highly accurate protein structure prediction with AlphaFold (Nature, 2021)
50. Abramson et al.: Accurate structure prediction of biomolecular interactions with AlphaFold 3 (Nature 630, 2024)
51. AlphaFold 3 inference pipeline and weights terms (Google DeepMind, 2026)
52. OpenFold3-preview repository (AlQuraishi Lab and OpenFold Consortium, 2025)
53. Boltz-2: Joint Structure and Binding Affinity Prediction (Boltz team, MIT Jameel Clinic and Recursion, 2025)
54. Accelerating scientific discovery with Co-Scientist (Google, arXiv, 2025)
55. Energy Department Launches Genesis Mission (U.S. DOE, 2025)
56. cuBLAS documentation: Floating Point Emulation (NVIDIA, 2026)
57. Unlocking Tensor Core Performance with Floating Point Emulation in cuBLAS (NVIDIA, 2025)
58. Environments (spack.yaml, spack.lock) (Spack documentation, 2026)
59. GPU and other device support (Apptainer documentation, 2026)
60. Jupyter on Open OnDemand (Open OnDemand, 2026)
61. Dask-Jobqueue documentation (Dask, 2026)
62. Deploying on Slurm (Ray documentation, 2026)