Vlasiator: Global Hybrid-Vlasov Simulator

• Vlasiator is a global hybrid-Vlasov plasma simulation code that models the solar wind–magnetosphere–ionosphere system by solving the Vlasov equation for ions and treating electrons as a fluid.
• It employs sparse velocity-space storage, semi-Lagrangian time integration, and adaptive mesh refinement to efficiently simulate six-dimensional phase-space dynamics.
• The code leverages advanced parallelism (MPI, OpenMP, GPU) and dynamic load balancing to bridge the gap between fluid and particle methods in large-scale magnetospheric studies.

Vlasiator is a global hybrid-Vlasov plasma simulation code for near-Earth space in which ions are treated kinetically through the Vlasov equation in six-dimensional phase space, while electrons are represented as a massless, charge-neutralizing fluid and the electromagnetic fields are evolved in a reduced Maxwell framework, commonly described through the Darwin or quasi-neutral approximation depending on the formulation adopted in a given study. It is designed to model the solar wind–magnetosphere–ionosphere system at ion-kinetic fidelity on supercomputer-scale architectures, combining a Cartesian spatial discretization, sparse velocity-space storage, semi-Lagrangian time integration, adaptive mesh refinement, and dynamic load balancing (Palmroth et al., 2018).

1. Concept and scope

Vlasiator is described as a “supercomputer-scale Vlasov theory–based model that incorporates the solar wind–magnetosphere–ionosphere system” and, in the global magnetospheric setting, solves a six-dimensional hybrid-Vlasov model on a spatially adaptive Cartesian grid where each spatial cell carries a sparse three-dimensional velocity-space block (Bouri et al., 2023, Kotipalo et al., 27 May 2025). In its standard formulation, ions are evolved through a phase-space distribution function, whereas electrons enter through a fluid closure; this separation reflects the stated design choice that modeling electrons kinetically would require a resolution orders of magnitude greater than ions, so ions are modeled kinetically and electrons as a fluid (Kotipalo, 2022).

The code has been used in both 2D–3V and 3D–3V configurations. The 2D–3V geometry appears in meridional-plane magnetospheric simulations with spatial resolution such as $\Delta x \approx 600\,\mathrm{km}$ and output cadence $\Delta t=1\,\mathrm{s}$ (Holmberg et al., 23 Sep 2025, Holmberg et al., 18 Jan 2026). Fully three-dimensional applications include global magnetospheric domains with adaptive refinement down to $\Delta x_{\min}=1000\,\mathrm{km}$ in regions such as the bow shock, magnetosheath, magnetopause, and tail current sheet (Grandin et al., 2023, Grandin et al., 2023).

A recurring motivation across the literature is that the model aims to bridge the gap between fluid descriptions and particle methods. The review literature states that it directly models the Vlasov equation on a grid instead of using the more common particle-in-cell approach, thereby avoiding statistical noise while retaining ion-kinetic dynamics at global scale (Palmroth et al., 2018, Battarbee et al., 2024).

2. Governing equations and hybrid closure

The ion dynamics are governed by the collisionless Vlasov equation. In the notation used repeatedly in the Vlasiator literature, the ion distribution function $f(\mathbf{x},\mathbf{v},t)$ satisfies

$$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$

with corresponding species-indexed variants when multiple ion species are considered [(Battarbee et al., 30 Apr 2025); (Palmroth et al., 2018); (Kempf et al., 2013)].

The electromagnetic evolution is closed by Faraday’s law and an Ampère relation without displacement current in the Darwin formulation,

$$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$

or, in some summaries, by the standard Vlasov–Maxwell form used to describe the underlying system solved by the code (Grandin et al., 2023, Holmberg et al., 23 Sep 2025, Bouri et al., 2023). The review literature further notes that in Vlasiator the “background” dipole and interplanetary magnetic field are held fixed and curl- and divergence-free, while the code advances the perturbation field $\delta\mathbf B$ (Palmroth et al., 2018).

Electrons are treated as a massless, charge-neutralizing fluid. Several closely related closures appear in the literature. A generalized Ohm’s law with Hall term and electron-pressure contribution is written as

$$\mathbf E = -\,\mathbf u_e\times\mathbf B +\frac{1}{e\,n_i}\,\mathbf J\times\mathbf B -\frac{1}{e\,n_i}\,\nabla p_e +\eta\,\mathbf J,$$

with $\mathbf u_e=\mathbf u_i-\mathbf J/(e\,n_i)$ in one formulation (Facskó, 2019). Other treatments give

$$\mathbf E = -\,\mathbf u_i\times\mathbf B +\frac{\mathbf J\times\mathbf B}{e n} -\frac{1}{e n}\nabla p_e,$$

with a polytropic closure $\Delta t=1\,\mathrm{s}$0 (Grandin et al., 2023), or emphasize the Hall-only reduction used in verification studies,

$\Delta t=1\,\mathrm{s}$1

contrasted with an ideal-MHD variant omitting the Hall term (Kempf et al., 2013). The review literature states that the Hall term is essential for recovering ion-kinetic resonances such as left- and right-hand mode splitting and Bernstein modes in the hybrid-Vlasov model (Kempf et al., 2013).

3. Phase-space discretization and numerical solution

Vlasiator discretizes the system on a Cartesian spatial grid and a Cartesian velocity grid carried by each spatial cell, but only a sparse subset of velocity space is stored. The review literature describes a user-set threshold $\Delta t=1\,\mathrm{s}$2 so that only velocity cells for which $\Delta t=1\,\mathrm{s}$3, plus one halo layer of neighbors in both $\Delta t=1\,\mathrm{s}$4 and $\Delta t=1\,\mathrm{s}$5, are stored and updated. The stored data are organized into small $\Delta t=1\,\mathrm{s}$6 velocity-cell sub-blocks to maximize data contiguity and vector-unit utilization (Palmroth et al., 2018). GPU-porting work reports that the global maximum velocities may be subdivided into blocks of size $\Delta t=1\,\mathrm{s}$7, historically $\Delta t=1\,\mathrm{s}$8 and later also $\Delta t=1\,\mathrm{s}$9, and that each spatial cell maintains a VelocityMesh and a VelocityBlockContainer for sparse storage and metadata (Battarbee et al., 2024).

Time integration is based on operator splitting. The six-dimensional Vlasov equation is split into spatial “translation” and velocity-space “acceleration” substeps using Strang splitting (Battarbee et al., 30 Apr 2025). The review literature describes a global leap-frog plus Strang-splitting cycle: an initial half-acceleration, followed by spatial translation $\Delta x_{\min}=1000\,\mathrm{km}$0, field update, and a second half-acceleration $\Delta x_{\min}=1000\,\mathrm{km}$1 (Palmroth et al., 2018). The semi-Lagrangian solver is based on SLICE-3D, which factorizes a full 3D remap into successive 1D remaps with conservative remapping; the magnetic-field-induced velocity-space rotation is further decomposed into orthogonal shears (Palmroth et al., 2018). The recent communication work specifies that Vlasiator uses a semi-Lagrangian SLICE-3D stencil of width 2 in each direction (Battarbee et al., 30 Apr 2025).

Adaptive mesh refinement is integral to the code’s spatial discretization strategy. Static targeted refinement uses the dccrg library and user-prescribed regions refined to level $\Delta x_{\min}=1000\,\mathrm{km}$2, subject to the rule that neighboring cells may differ by at most one refinement level (Kotipalo, 2022). Adaptive refinement based on simulation data employs dimensionless gradient indicators. For a set of variables including ion density $\Delta x_{\min}=1000\,\mathrm{km}$3, perturbed field energy density $\Delta x_{\min}=1000\,\mathrm{km}$4, momentum density $\Delta x_{\min}=1000\,\mathrm{km}$5, and perturbed magnetic field $\Delta x_{\min}=1000\,\mathrm{km}$6, the refinement indicator is defined as

$\Delta x_{\min}=1000\,\mathrm{km}$7

with refinement and coarsening controlled by thresholds $\Delta x_{\min}=1000\,\mathrm{km}$8 and $\Delta x_{\min}=1000\,\mathrm{km}$9 (Kotipalo, 2022). In reported tests, AMR with two levels of refinement and approximate thresholds $f(\mathbf{x},\mathbf{v},t)$0 and $f(\mathbf{x},\mathbf{v},t)$1 concentrated resolution near bow-shock flanks, reduced CPU time per global timestep by $f(\mathbf{x},\mathbf{v},t)$2 compared to static refinement, and used $f(\mathbf{x},\mathbf{v},t)$3 the coarse-mesh memory rather than the $f(\mathbf{x},\mathbf{v},t)$4 of static refinement (Kotipalo, 2022).

4. Parallel architecture, communication, and load balancing

Vlasiator is built for distributed-memory high-performance computing. Spatial domain decomposition is handled through the DCCRG library, with dynamic rebalancing through Zoltan (Battarbee et al., 30 Apr 2025, Palmroth et al., 2018). The literature describes three nested levels of parallelism: MPI domain decomposition in physical space, OpenMP threading within each MPI domain, and SIMD vectorization inside the 1D remaps operating on the velocity blocks (Palmroth et al., 2018). GPU-porting work extends this to heterogeneous CPU/GPU execution through a unified CUDA/HIP codebase, Unified Memory for core phase-space structures, GPU streams controlled by OpenMP threads, and on-device sparse hash maps for velocity-mesh operations (Battarbee et al., 2024).

Load balancing is a central issue because the number of active velocity-space blocks and the spatial refinement level vary by orders of magnitude across the domain (Kotipalo et al., 27 May 2025). The 2025 load-balancing study formulates the problem as weighted graph or hypergraph partitioning. In the graph model, each vertex corresponds to a spatial cell with computational weight equal to the number of active velocity-space blocks, and each edge carries a communication weight based on the same block count when neighboring cells lie on different ranks (Kotipalo et al., 27 May 2025). Four heuristic families are examined in Vlasiator: parallel hypergraph partitioner (PHG), recursive coordinate bisection (RCB), recursive inertial bisection (RIB), and Hilbert space-filling curve (HSFC) partitioning (Kotipalo et al., 27 May 2025).

The reported large-scale tests on EuroHPC LUMI compare a sample run with 285 972 spatial cells and $f(\mathbf{x},\mathbf{v},t)$5 total velocity blocks against a large run with 11 940 800 spatial cells and $f(\mathbf{x},\mathbf{v},t)$6 blocks (Kotipalo et al., 27 May 2025). In those tests, HSFC with the Beta curve provided the best overall performance. The measured post-repartition load-imbalance $f(\mathbf{x},\mathbf{v},t)$7 on the large case was $f(\mathbf{x},\mathbf{v},t)$8 for RCB, $f(\mathbf{x},\mathbf{v},t)$9 for RIB, $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$0 for graph and hypergraph PHG variants, and $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$1 for HSFC (Beta) (Kotipalo et al., 27 May 2025). On the large case, HSFC (Beta) outperformed RIB by $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$2 in spatial propagation and $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$3 overall, while among Hilbert variants the Beta curve outperformed the default Octree ordering by $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$4 in spatial propagation and $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$5 overall (Kotipalo et al., 27 May 2025). Earlier work on BlueGene/P had already identified partitioning overhead as a substantial component of total runtime and found that graph partitioners could achieve nearly the same balancing quality as PHG with lower preprocessing cost (Turk et al., 2011).

A separate communication optimization introduced “Ghost Translation,” a coalesced halo-exchange strategy for the translation step of the 6D Vlasov solve (Battarbee et al., 30 Apr 2025). In the classical solver, six separate MPI exchange steps per translation are used; Ghost Translation replaces these with a single halo exchange that gathers all needed ghost-cell data up front (Battarbee et al., 30 Apr 2025). The cost model is written as

$$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$6

where $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$7 is per-message latency, $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$8 inverse bandwidth, $$\frac{\partial f}{\partial t} + \mathbf v\cdot\nabla_{\mathbf x}f + \frac{q}{m}\bigl(\mathbf E+\mathbf v\times\mathbf B\bigr)\cdot\nabla_{\mathbf v}f =0,$$9 the size of one face halo, and $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$0 classically (Battarbee et al., 30 Apr 2025). On up to 8 000 ranks on the Mahti cluster, the minimal-ghost variant GT1 cut MPI latencies by a factor of 3–5 in weak scaling and, in a realistic 8 000-rank translate-only test, changed compute time from $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$1 to $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$2, MPI communication from $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$3 to $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$4, and memory footprint by $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$5 relative to the classical method; GT3 increased memory by $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$6 and became much heavier when AMR enlarged active regions (Battarbee et al., 30 Apr 2025).

5. Verification and physical interpretation

Verification of Vlasiator has relied in part on wave-dispersion studies against WHAMP. In one-dimensional low-$$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$7 plasma tests with $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$8, $$\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \nabla\times\mathbf B=\mu_0\mathbf J,\qquad \nabla\cdot\mathbf B=0,$$9, $\delta\mathbf B$0, domain length $\delta\mathbf B$1 with 10 000 cells, velocity-space spacing $\delta\mathbf B$2, and time step $\delta\mathbf B$3, the hybrid-Vlasov dispersion curves agreed well with WHAMP in the low-frequency, long-wavelength regime relevant for global magnetospheric simulations (Kempf et al., 2013). The reported results include agreement of left- and right-hand polarized modes and ion Bernstein modes, while also showing that omission of the Hall term reduces the response to non-dispersive MHD Alfvén behavior (Kempf et al., 2013).

The broader review literature lists additional verification and validation elements, including comparison of low-$\delta\mathbf B$4 wave dispersion with WHAMP, ion/ion right-hand resonant beam instability growth rates, comparison to hybrid-PIC in one-dimensional setups, and 1D collisionless shock studies (Palmroth et al., 2018). It also describes validations against spacecraft and benchmark models, including foreshock ion beams and $\delta\mathbf B$5 ULF waves, mirror-mode structures downstream of a quasi-perpendicular shock, dayside flux transfer event generation, and tail reconnection with plasmoid ejection (Palmroth et al., 2018). This suggests that Vlasiator has been developed not only as a numerical platform but also as a physics model whose credibility rests on dispersion, shock, and mesoscale magnetospheric tests.

At the same time, the literature records critical interpretation issues. In commentary on simulations of foreshock transients, it is argued that events labeled as spontaneous hot flow anomalies lacked expected strong heating, characteristic magnetic turbulence, pronounced solar-wind flow deflection and deceleration, and the “shoulders” of density and field enhancement at cavity edges; those simulated structures are instead interpreted as proto-SHFAs (Facskó, 2019). The stated limitations include insufficient spatial and velocity-space resolution to capture thin boundary layers, possible underprediction of ion heating in under-resolved regions, and global setups in which structures may advect out of the simulation volume before fully nonlinear development (Facskó, 2019).

6. Scientific applications and current extensions

Vlasiator has been used as a physics engine for diverse downstream analyses. In soft X-ray modeling, a 3D–3V global simulation with base cell $\delta\mathbf B$6, refinement to $\delta\mathbf B$7, and velocity resolution $\delta\mathbf B$8 was post-processed to compute solar-wind charge-exchange emissivity and line-of-sight integrated images, yielding signatures of mirror-mode structures in the magnetosheath and flux transfer events at the magnetopause (Grandin et al., 2023). In auroral precipitation studies, the first 3D hybrid-Vlasov global simulation of auroral proton precipitation reproduced cusp bursts associated with flux transfer events and nightside injections associated with fast Earthward plasma flows in the near-Earth magnetotail, with simulated precipitating fluxes reported to be in good agreement with DMSP observations under comparable driving (Grandin et al., 2023).

The code has also become a source of high-dimensional datasets for topology and machine learning. One magnetic-topology study analyzed Vlasiator magnetotail outputs on a $\delta\mathbf B$9 grid with 1 000 km spacing and 1 s cadence, detecting 2D and 3D magnetic nulls and constructing spatial and spatiotemporal graphs of null connectivity (Bouri et al., 2023). Subsequent graph-based surrogate models used Vlasiator outputs to emulate near-Earth space dynamics. One 2025 study trained deterministic and probabilistic graph neural models on 2D–3V Vlasiator data and reported next-step inference times of about $$\mathbf E = -\,\mathbf u_e\times\mathbf B +\frac{1}{e\,n_i}\,\mathbf J\times\mathbf B -\frac{1}{e\,n_i}\,\nabla p_e +\eta\,\mathbf J,$$0 for deterministic forecasts and about $$\mathbf E = -\,\mathbf u_e\times\mathbf B +\frac{1}{e\,n_i}\,\mathbf J\times\mathbf B -\frac{1}{e\,n_i}\,\nabla p_e +\eta\,\mathbf J,$$1 for a probabilistic five-member ensemble on one GPU, compared with 4–5 min per simulated second for the original model on 50 CPUs (Holmberg et al., 23 Sep 2025). A later study on a 670k-cell grid reported more than two orders of magnitude speedup relative to runs on 100 CPUs, while incorporating a magnetic-divergence penalty and CRPS fine-tuning for uncertainty-aware forecasting (Holmberg et al., 18 Jan 2026). These studies do not alter Vlasiator’s numerical core, but they indicate that its outputs are being treated as a reference dynamical system for real-time emulation.

Recent methodological development has extended beyond performance engineering into coupled geospace modeling. An ionosphere-coupled Vlasiator configuration maps field-aligned current density, density, temperature, and precipitating-electron energy flux from a coupling radius near $$\mathbf E = -\,\mathbf u_e\times\mathbf B +\frac{1}{e\,n_i}\,\mathbf J\times\mathbf B -\frac{1}{e\,n_i}\,\nabla p_e +\eta\,\mathbf J,$$2 to a thin-shell ionosphere at $$\mathbf E = -\,\mathbf u_e\times\mathbf B +\frac{1}{e\,n_i}\,\mathbf J\times\mathbf B -\frac{1}{e\,n_i}\,\nabla p_e +\eta\,\mathbf J,$$3, solves

$$\mathbf E = -\,\mathbf u_e\times\mathbf B +\frac{1}{e\,n_i}\,\mathbf J\times\mathbf B -\frac{1}{e\,n_i}\,\nabla p_e +\eta\,\mathbf J,$$4

and feeds the resulting ionospheric electric field back to the magnetospheric boundary (Horaites et al., 2 Feb 2026). In that framework, FTEs at the magnetopause were related to geoelectric-field structures near their magnetic footpoints, with 3D magnetic null points separating individual flux ropes in the simulation (Horaites et al., 2 Feb 2026). A plausible implication is that Vlasiator is evolving from a global hybrid-Vlasov magnetospheric solver into a more explicitly coupled geospace platform, while retaining the same foundational emphasis on ion kinetics, sparse phase-space representation, and supercomputer-scale parallelism.