Massive stochastic simulation of cosmic rays propagation in the heliosphere: The COSMICA code

The accurate modeling of galactic cosmic ray (GCR) propagation in the heliosphere requires solving the Parker Transport Equation (PTE), a multidimensional nonlinear equation that cannot be addressed analytically without strong approximations. In recent decades, stochastic differential equation (SDE)–Monte Carlo methods have emerged as a powerful numerical strategy for this problem, thanks to their numerical stability, relatively low memory requirements, and intrinsic parallelism. The increasing availability of general-purpose Graphics Processing Units (GPUs) has further revolutionized this approach by enabling massive parallelization of particle trajectories at relatively low cost. In this work, we introduce COSMICA (COde for a Speedy Montecarlo Involving Cuda Architecture), a new open-source multi-GPU code written in CUDA/C++ for the three-dimensional solution of the PTE. COSMICA has been specifically designed to optimize GPU resource usage and scalability, with strategies including memory hierarchy exploitation, register-conscious kernel design, warp-aware scheduling, and parameter reordering for multi-GPU execution. Benchmark results demonstrate that COSMICA reduces runtimes from weeks to hours for large-scale simulations. These optimizations make COSMICA a versatile tool for systematic studies of cosmic-ray modulation and parameter exploration, thereby expanding the feasibility of investigations that were previously computationally prohibitive. The present article constitutes the first part of a two-paper series, focusing on code design and computational performance; a companion paper will present its validation against benchmark models.