Immersed boundary methods (IBMs) facilitate the simulation of flows around stationary, moving, and deforming bodies on Cartesian grids. However, extending these simulations to the large grid sizes required for realistic flow problems remains a significant computational challenge. In this work, we present the implementation and acceleration of ViCar3D, a sharp-interface immersed boundary solver, on graphical processing units (GPUs). We utilize OpenACC, CUDA Fortran and MPI to reprogram ViCar3D, a sharp-interface immersed boundary solver, on multi-GPU architectures. Verification and scalability studies are performed for two benchmark cases: two-dimensional flow past a circular cylinder and direct numerical simulation (DNS) of flow past a finite rectangular wing. For the latter, we observe an approximately 20X speedup (node-to-node comparison) relative to the CPU-based implementation. The GPU-accelerated solver is capable of simulating complex 3D flows with up to 200 million mesh points on a single node equipped with four GPUs. Strong and weak scaling tests demonstrate maximum scaling efficiencies of 92% and 93%, respectively, on multi-GPU systems. We further test the code to simulate fluid flow past complex-shaped single-body and multi-body cases.

We present a scalable incompressible Navier–Stokes solver for three-dimensional Cartesian grids with non-uniform spacing. A direct tensor-product–Thomas method solves the constant-coefficient Poisson and Helmholtz equations arising from pressure projection and implicit diffusion, without approximate factorization. On uniform grids, the method recovers the classical eigenfunction-expansion method evaluated with fast Fourier transforms (FFTs). On stretched grids, diagonal scaling symmetrizes the one-dimensional Laplace operators, and the resulting numerical eigenbasis transforms are evaluated as general matrix–matrix multiplications (GEMMs). FFTs and GEMMs can be selected independently in each diagonalized direction while retaining the pencil decomposition, collective transposes, and tridiagonal machinery of an established FFT-based solver. The elliptic solver is verified to round-off accuracy, and the complete flow solver is validated against benchmark flows. Against geometric multigrid and block cyclic reduction with FFT diagonalization, the present method achieves the lowest time-to-solution among the tested approaches; at fixed grid dimensions, its cost is insensitive to grid stretching. CPU and multi-GPU tests show that GEMM-rich variants attain higher strong-scaling efficiency by better amortizing communication. On a single GPU, the fully GEMM-based variant increases the Poisson cost by 2.8× but the complete Navier–Stokes step cost by only 1.8×. Weak scaling exposes the trade-off: dense-transform costs grow with the global transform dimension, whereas non-uniform meshes can reduce the required number of grid points. The resulting open-source solver, , extends an FFT-based Navier–Stokes solver to large-scale simulations on grids stretched in multiple directions.

Computational fluid dynamics and fluid-structure interaction simulations involving moving and deforming bodies is extremely hard. In this work, we present a graphical processing unit (GPU) optimized implementation of the sharp-interface immersed boundary method. The method allows performing simulation around complex stationary as well as moving bodies on a Cartesian grid. We base our implementation on the ViCar3D framework and make use of OpenACC, CUDA, NCCL and MPI. We test the implementation across grid sizes ranging from O(10million) to O(1billion) points and achieved a 20X speedup compared to existing CPU implementation. We next present our multi-GPU implementation by utilizing CUDA streams and NCCL communicators. This enables us to obtain a >90

We explicitly demonstrate the universality of critical dynamics through unprecedented large-scale Graphics Processing Units (GPU)-based simulations of two out-of-equilibrium processes, comparing the behavior of spin-1/2 Ising and spin-1 Blume-Capel models on a square lattice. In the first protocol, a completely disordered system is instantaneously brought into contact with a thermal bath at the critical temperature, allowing it to evolve until the coherence length exceeds 10^{3} lattice spacings. Finite-size effects are negligible due to the mesoscopic scale of the lattice sizes studied, with linear dimensions up to L=2^{22} and 2^{19} for the Ising and Blume-Capel models, respectively. Our numerical data, and the subsequent analysis, demonstrate a strong dynamic universality between the two models and provide the most precise estimate to date of the dynamic critical exponent for this universality class, z=2.1676(1). In the second protocol, we corroborate the role of the universal ratio of dynamic and static length scales in achieving an exponential acceleration in the approach to equilibrium just above the critical temperature, through a time-dependent variation of the thermal bath temperature. The results presented in this work leverage our Compute Unified Device Architecture (CUDA)-based numerical code, breaking the world record for the simulation speed of the Ising model.

We present high-performance implementations of the two-dimensional Ising and Blume-Capel models for largescale, multi-GPU simulations. Our approach takes full advantage of the NVIDIA GB200 NVL72 system, which features up to 72 GPUs interconnected via high-bandwidth NVLink, enabling direct GPU-to-GPU memory access across multiple nodes. By utilizing Fabric Memory and an optimized Monte Carlo kernel for the Ising model, our implementation supports simulations of systems with linear sizes up to L = 2(23), corresponding to approximately 70 trillion spins. This allows for a peak processing rate of nearly 1.15 x 10(5) lattice updates per nanosecond- setting a new performance benchmark for Ising model simulations. Additionally, we introduce a custom protocol for computing correlation functions, which strikes an optimal balance between computational efficiency and statistical accuracy. This protocol enables large-scale simulations without incurring prohibitive runtime costs. Benchmark results show near-perfect strong and weak scaling up to 64 GPUs, demonstrating the effectiveness of our approach for large-scale statistical physics simulations. Program summary Program title: cuIsing (optimized) CPC Library link to program files: https://doi.org/10.17632/ppkwwmcpwg.1 Licensing provisions: MIT license Programming languages: CUDA C Nature of problem: Comparative studies of the critical dynamics of the Ising and Blume-Capel models are essential for gaining deeper insights into phase transitions, enhancing computational methods, and developing more accurate models for complex physical systems. To minimize finite-size effects and optimize the statistical quality of simulations, large-scale simulations over extended time scales are necessary. To support this, we provide two high-performance codes capable of running simulations with up to 70 trillion spins. Solution method: We present updated versions of our multi-GPU code for Monte Carlo simulations, implementing both the Ising and Blume-Capel models. These codes take full advantage of multi-node NVLink systems, such as the NVIDIA GB200 NVL72, enabling scaling across GPUs connected across different nodes within the same NVLink domain. Communication between GPUs is handled seamlessly via Fabric Memory-a novel memory allocation technique that facilitates direct memory access between GPUs within the same domain, eliminating the need for explicit data transfers. By employing highly optimized CUDA kernels for the Metropolis algorithm and a custom protocol that reduces the computational overhead of the correlation function, our implementation achieves the highest recorded performance to date.

COVID-19 has shown the importance of having a fast response against pandemics. Finding a novel drug is a very long and complex procedure, and it is possible to accelerate the preliminary phases by using computer simulations. In particular, virtual screening is an in-silico phase that is needed to filter a large set of possible drug candidates to a manageable number. This paper presents the implementations and a comparative analysis of two GPU-optimized implementations of a virtual screening algorithm targeting novel GPU architectures. The first adopts a traditional approach that spreads the computation required to evaluate a single molecule across the entire GPU. The second uses a batched approach that exploits the parallel architecture of the GPU to evaluate more molecules in parallel, without considering the latency to process a single molecule. The paper describes the advantages and disadvantages of the proposed solutions, highlighting implementation details that impact the performance. Experimental results highlight the different performance of the two methods on several target molecule databases while running on NVIDIA A100 GPUs. The two implementations have a strong dependency with respect to the data to be processed. For both cases, the performance is improving while reducing the dimension of the target molecules (number of atoms and rotatable bonds). The two methods demonstrated a different behavior with respect to the size of the molecule database to be screened. While the latency one reaches sooner (with fewer molecules) the performance plateau in terms of throughput, the batched one requires a larger set of molecules. However, the performances after the initial transient period are much higher (up to 5x speed-up). Finally, to check the efficiency of both implementations we deeply analyzed their workload characteristics using the instruction roof-line methodology.

This chapter shows some applications of the Fast Fourier Transform and how to utilize the CUFFT library.

Diffusion tensor imaging (DTI) is a non-invasive magnetic resonance imaging (MRI) modality used to map white matter fiber tracts for a variety of clinical applications; one of which is aiding preoperative assessments for tumor patients. DTI requires numerical computations on multiple diffusion weighted images to calculate diffusion tensors at each voxel and probabilistic tracking1 to construct fiber tracts, or tractography. Greater accuracy in tractography is possible with larger, more advanced imaging and reconstruction algorithms. However, larger scans and advanced reconstruction is often computationally intensive. The post-processing pipeline involves significant computational resources and time and requires up to 40 minutes of computation time on state-of-the-art hardware. Parallel GPU computations can improve time for the resource-intensive tractography. A collaborative team from DIPY, NVIDIA, and UCSF recently developed a tool, GPUStreamlines, for GPU-enabled tractography2 which has been expanded to support the constant solid angle (CSA) reconstruction algorithm3. This GPU-enabled tractography was applied to MRIs of brains with and without presence of lesions, with substantial increases in processing speed. We demonstrate that CSA GPU-enabled tractography in normal controls and patients are comparable to the existing gold standard tractography currently in place at UCSF.

Virtual screening is an early stage of the drug discovery process that selects the most promising candidates. In the urgent computing scenario it is critical to find a solution in a short time frame. In this paper, we focus on a real-world virtual screening application to evaluate out-of-kernel optimizations, that consider input and architecture features to improve the computation efficiency on GPU. Experiment results on a modern supercomputer node show that we can almost double the performance. Moreover, we implemented the optimization using SYCL and it provides a consistent benefit with the CUDA optimization. A virtual screening campaign can use this gain in performance to increase the number of evaluated candidates, improving the probability of finding a drug.

This paper is associated with a video winner of a 2021 American Physical Society's Division of Fluid Dynamics (DFD) Gallery of Fluid Motion Award for work presented at the DFD Gallery of Fluid Motion. The original video is available online at the Gallery of Fluid Motion, https://doi.org/10.1103/APS.DFD.2021.GFM.V0053.

This paper presents a performance analysis of pencil domain decomposition methodologies for three-dimensional Computational Fluid Dynamics (CFD) codes for turbulence simulations, on several large GPU-accelerated clusters. The performance was assessed for the numerical solution of the Navier-Stokes equations in two codes which require the calculation of Fast-Fourier Transforms (FFT): a tri-periodic pseudo-spectral solver for isotropic turbulence, and a finite-difference solver for canonical turbulent flows, where the FFTs are used in its Poisson solver. Both codes use a newly developed transpose library that automatically determines the optimal domain decomposition and communication backend on each system. We compared the performance across systems with very different node topologies and available network bandwidth, to show how these characteristics impact decomposition selection for best performance. Additionally, we assessed the performance of several communication libraries available on these systems, such as Open-MPI, IBM Spectrum MPI, Cray MPI, the NVIDIA Collective Communication Library (NCCL), and NVSHMEM. Our results show that the optimal combination of communication backend and domain decomposition is highly system-dependent, and that the adaptive decomposition library is key in ensuring efficient resource usage with minimal user effort.

The reliability of cardiovascular computational models depends on the accurate solution of the hemodynamics, the realistic characterization of the hyperelastic and electric properties of the tissues along with the correct description of their interaction. The resulting fluid-structure-electrophysiology interaction (FSEI) thus requires an immense computational power, usually available in large supercomputing centers, and requires long time to obtain results even if multi-CPU processors are used (MPI acceleration). In recent years, graphics processing units (GPUs) have emerged as a convenient platform for high performance computing, as they allow for considerable reductions of the time-to-solution. This approach is particularly appealing if the tool has to support medical decisions that require solutions within reduced times and possibly obtained by local computational resources. Accordingly, our multi-physics solver has been ported to GPU architectures using CUDA Fortran to tackle fast and accurate hemodynamics simulations of the human heart without resorting to large-scale supercomputers. This work describes the use of CUDA to accelerate the FSEI on heterogeneous clusters, where both the CPUs and GPUs are used in synergistically with minor modifications of the original source code. The resulting GPU accelerated code solves a single heartbeat within a few hours (from three to ten depending on the grid resolution) running on premises computing facility made of few GPU cards, which can be easily installed in a medical laboratory or in a hospital, thus opening towards a systematic computational fluid dynamics (CFD) aided diagnostic.

We study the statistics of passive scalars at $Pr=1$ , for turbulent flow within a smooth straight pipe of circular cross section up to $Re_{\tau } \approx 6000$ using direct numerical simulation (DNS) of the Navier–Stokes equations. While featuring a general organisation similar to the axial velocity field, passive scalar fields show additional energy at small wavenumbers, resulting in a higher degree of mixing and in a $k^{-4/3}$ spectral inertial range. The DNS results highlight logarithmic growth of the inner-scaled bulk and mean centreline scalar values with the friction Reynolds number, implying an estimated scalar von Kármán constant $k_{\theta } \approx 0.459$ , which also nicely fits the mean scalar profiles. The DNS data are used to synthesise a modified form of the classical predictive formula of Kader & Yaglom (Intl J. Heat Mass Transfer, vol. 15 (12), 1972, pp. 2329–2351), which points to some shortcomings of the original formulation. Universality of the mean core scalar profile in defect form is recovered, with very nearly parabolic shape. Logarithmic growth of the buffer-layer peak of the scalar variance is found in the Reynolds number range under scrutiny, which well conforms with Townsend's attached-eddy hypothesis, whose validity is also supported by the spectral maps. The behaviour of the turbulent Prandtl number shows good universality in the outer wall layer, with values $Pr_t \approx 0.84$ , as also found in previous studies, but closer to unity near the wall, where existing correlations do not reproduce the observed trends.