Nyx: Adaptive N-body & Hydro Code

Updated 21 December 2025

Nyx is a massively parallel simulation framework that integrates a symplectic N-body solver with a second-order unsplit Godunov hydrodynamic scheme on adaptive mesh grids.
It employs a block-structured AMR framework with advanced multigrid methods to accurately model cosmological structure formation and gravitational dynamics.
The Nyx-RT extension couples adaptive, photon-conserving ray tracing with gas and dark matter dynamics, enabling efficient radiative transfer on exascale platforms.

Nyx is a massively parallel block-structured adaptive mesh refinement (AMR) code for large-scale cosmological N-body and hydrodynamics simulations. It integrates a symplectic particle-mesh N-body solver for dark matter, a second-order unsplit Godunov Eulerian scheme for gas dynamics, and robust Poisson gravity on dynamic grids. Its architecture is built for exascale, leveraging hybrid MPI+OpenMP/GPU parallelization and the AMReX framework. Recent extensions, notably Nyx-RT, incorporate fully adaptive, photon-conserving ray tracing for radiative transfer, coupling radiation with gravitational and gas dynamics efficiently at unprecedented resolution and scale (Almgren et al., 2013, Marshak et al., 13 Dec 2025).

1. Mathematical and Algorithmic Foundations

Nyx solves cosmological structure formation by evolving dark matter and baryonic gas in an expanding universe. The comoving scale factor $a(t)$ is advanced using the Friedmann equation:

$\frac{d(\ln a)}{dt} = H_0 \sqrt{\Omega_0 a^{-3} + \Omega_\Lambda}$

Time integration of $a(t)$ uses a second-order Runge–Kutta method, maintaining $|\Delta a/a| < 10^{-8}$ per step.

Gas dynamics are formulated in comoving coordinates, advancing the conservation laws for baryon density $\rho_b$ and peculiar velocity $U$ :

$\partial_t \rho_b = -\frac{1}{a} \nabla \cdot (\rho_b U) \ \partial_t (a \rho_b U) = -\nabla \cdot (\rho_b U U) - \nabla p + \rho_b g$

The dual-energy formalism evolves both total energy $E$ and internal energy $e$ for accuracy in hypersonic and low-Mach regimes, with:

$\partial_t (a^2 \rho_b E) = - a\nabla\cdot (\rho_b U E + p U) + a \rho_b U \cdot g + a \Lambda_{HC}$

The N-body module replaces the collisionless Boltzmann equation with $N$ particles, each evolved using a symplectic kick–drift–kick algorithm that naturally incorporates expansion:

// Pseudocode: particle update
(a u_i)^{n+1/2} = (a u_i)^n + (Δt/2) g^n(x_i)
x_i^{n+1} = x_i^n + (Δt / a_{n+1/2}) u_i^{n+1/2}
(a u_i)^{n+1} = (a u_i)^{n+1/2} + (Δt/2) g^{n+1}(x_i^{n+1})

Self-gravity is computed via

$\nabla^2 \phi = \frac{4\pi G}{a} (\rho_b + \rho_{dm} - \bar\rho )$

with periodic boundaries and densities deposited using the cloud-in-cell (CIC) method (Almgren et al., 2013).

2. Numerical Methods and AMR Framework

Nyx’s block-structured AMR allows spatial refinement by factors of two, utilizing the Berger–Rigoutsos clustering algorithm for grid generation, and supporting multiple subcycling policies for temporal refinement. Algorithms for mesh and multilevel synchronization include:

Proper nesting: each fine patch contained within coarser regions with a buffer border.
Recursive multilevel advance: gravity is solved on all levels, hydrodynamics advanced, particle update performed, and synchronization executed via refluxing and elliptic corrections.

The Eulerian hydro solver is an unsplit Godunov method, employing either a piecewise linear method (PLM) or enhanced piecewise parabolic method (PPM) for reconstruction, and an iterative two-shock (Colella–Glaz) Riemann solver for robustness in shocks.

Mass assignment and force interpolation in the particle-mesh N-body algorithm use CIC across the AMR hierarchy. Ghost and virtual particles ensure correct mass mapping at coarse–fine boundaries. The Poisson equation is discretized with a 7-point stencil and solved by geometric multigrid using red-black Gauss–Seidel smoothing.

3. Parallel Implementation and Workflow

Nyx is written in Fortran90 and C++ and leverages hybrid MPI+OpenMP parallelism (BoxLib/AMReX). Key aspects include:

Domain decomposition: per-level “MultiFab” grids assigned via Knapsack heuristics or Morton-order (Z-curve) for load balance.
Fine-grained data partitioning: Each MPI rank manages its particles and local mesh blocks, with no global particle structure.
Communication: Ghost-cell halo exchange for hydro and gravity, global redistributions only on regridding or grid boundary crossing. Multigrid boundary data is communicated for solution consistency.
I/O: Output is stored as self-describing HDF5 “plotfile” and “checkpoint” directories, with coordinated job restart and optional byte-swap or precision conversion. Built-in support exists for amrvis, VisIt, and yt visualization tools (Almgren et al., 2013).

4. Radiative Transfer: Nyx-RT Adaptive Ray Tracing

The Nyx-RT extension provides operator-split coupling of radiative transfer (RT) with gravitational and gas dynamics. Adaptive ray tracing is implemented as a sequence of AMReX particle operations, where each “ray” is a pseudo-particle carrying photon luminosity, travel distance, optical depth, and direction (parameterized on a HEALPix hierarchy). The algorithm includes:

Initializing rays as level-0 HEALPix bundles at each source.
Marching rays cell-by-cell using Amanatides–Woo voxel traversal.
Depositing photoionization and photoheating rates locally, with photon-conserving flux calculation:

$\Gamma = \int_{\nu_0}^{\infty} \frac{L_\nu e^{-\tau_\nu}}{h\nu} \frac{1-e^{-\Delta\tau_\nu}}{n_{\text{HI}} V_{\text{cell}}} d\nu$

Splitting rays adaptively as determined by local solid angle coverage; filtering rays based on optical depth, travel distance, or pixel depth.

Geometric overlap corrections mitigate numerical artifacts from ray–pixel mismatch, reassigning “lost” photons to neighboring cells with a conservative thresholding procedure for numerical stability (Marshak et al., 13 Dec 2025).

Source merging introduces dynamic consolidation of ionizing sources whose H II regions overlap, using a radius-detection and union-find graph algorithm. This reduces the number of rays by $\sim60\%$ in late reionization, halving the total wall-time for 1024³ runs with only minimal effects ( $|\Delta z|{\lesssim}0.2$ ) on reionization history.

5. Validation and Benchmarking

Nyx and Nyx-RT have undergone extensive code validation:

N-body particle dynamics: Two-particle orbital tests show radius drift $<1.1\%$ (coarse) and $<0.17\%$ (fine), with multi-level accuracy matching single-level within $2\times10^{-4}$ .
Cosmological N-body: Two-point correlation $\xi(r)$ , halo mass functions, and NFW profiles match Gadget-2 and other AMR codes, with AMR runs (256³+2 levels) approximating uniform 1024³ grid results.
Hydrodynamical tests: Adiabatic Santa Barbara cluster simulations recover $M_{200}$ and $r_{200}$ , and reproduce radial profiles consistent with Enzo/ART and code-averaged benchmarks, including characteristic entropy core flattening.
Radiative transfer: Strömgren sphere solutions converge to analytic expectations; cosmological “Test 4” matches established results (Arc). Figures 4–5 in (Marshak et al., 13 Dec 2025) illustrate agreement for ionized regions and artifact suppression.

6. Performance and Scalability

Nyx is architected for exascale platforms, with demonstrated scaling to 49,152 CPU cores in MPI+OpenMP and to 4,096 GPUs with Nyx-RT. Scaling results show:

Weak scaling with $<50\%$ time/step increase over $1,000\times$ problem size expansion.
At 4096³ grid and particle resolution, full RT+hydrodynamics simulations are tractable ( $\mathcal{O}(10^4)$ GPU-hours for $z\sim199\to5$ ), with the ray-tracing pipeline GPU-dominated for memory bandwidth and throughput.
Source merging reduces total ray count by more than threefold, dramatically accelerating late-stage reionization calculations (Almgren et al., 2013, Marshak et al., 13 Dec 2025).

7. Applications, Impact, and Future Directions

Nyx and Nyx-RT are targeted at precision studies of cosmological structure formation, the epoch of reionization, and IGM properties. Demonstrated applications include:

Large-volume cosmological N-body/hydro with full radiative feedback at 4096³ scale.
Accurate Ly $\alpha$ forest flux modeling with convergence at the 10% level up to $k\sim1\,h\,\mathrm{Mpc}^{-1}$ .
Flexible AMR capabilities for multi-scale galaxy formation and cluster studies.

Ongoing and future development areas include improved photoheating models, helium radiative transfer, device-level GPU optimizations, ray merging, and full extension of adaptive RT to AMR (block-refined) grids (Marshak et al., 13 Dec 2025).

Summary Table: Key Capabilities of Nyx and Nyx-RT

Capability	Implementation	Benchmark Highlights
N-body solver	Particle-mesh, CIC, symplectic KDK	$\xi(r)$ , mass function, NFW
Hydrodynamics	Unsplit Godunov, PPM/PLM, dual-energy	Santa Barbara cluster, shocks
Gravity	Geometric multigrid Poisson	Scalable to > $10^4$ cores/GPUs
AMR	Block-structured, flexible subcycling	Multi-level accuracy, load balance
Radiative transfer (RT)	Adaptive ray tracing, source merging	Strömgren, full cosmological RT