MUSIC 4 GPU (MUSIG)

MUSIC is a 3+1D relativistic second-order viscous hydrodynamics for heavy ion collisions written in C++: https://github.com/MUSIC-fluid/MUSIC. A user manual for MUSIC is available here: https://webhome.phy.duke.edu/~jp401/old_music_manual/

---

GPU ports — scope and status (read first)

This repository adds GPU acceleration (CUDA and Metal back-ends) on top of MUSIC. Before using it, understand exactly what it is and is not:

• No new physics. The GPU code is a direct port of the existing CPU algorithms in MUSIC release 3.1 — identical equations and numerical scheme (KT flux, Newton reconstruction, second-order Israel–Stewart viscous evolution). Nothing physical was added, removed, or changed.
• Authorship. The GPU ports, benchmarks, and GPU documentation were produced with Claude Opus 4.7.

✅ VALIDATED AGAINST THE CPU BUILD — A LARGE-SCALE CAMPAIGN COMPARISON IS STILL MISSING

The GPU paths have been compared with the CPU build (MUSIC_FORCE_CPU=1, the same binary, double precision) at every level from the evolution to final-state hadrons:

• Numerics: eps_max within 1e-4 (EOS 0, Gubser) and 6.4e-4 (hotQCD, EOS 91); one O+O event (3D Glauber) with e, vx, vy, vz within 1e-3.
• A 0–10% Au+Au jet production (X-SCAPE: 3D MC-Glauber strings, a viscous background and a jet leg with liquefier droplets, iSS on the freeze-out surfaces; 6 jet events on 3 backgrounds, CUDA on a GB10), see docs/MUSIC_CPU_vs_GPU.md: identical freeze-out times in all 9 hydro runs; energy-weighted |Δe|/e = 2×10⁻³ and |ΔT|/T = 3×10⁻⁴ above freeze-out; momentum anisotropy, ⟨pT⟩, v₂, v₃ and the jet wake unchanged within errors. Systematic: on the GPU 0.3% less energy above freeze-out, a 0.19% smaller freeze-out volume and 0.20 ± 0.02% fewer charged hadrons.
• Metal (M3 Max): the same production test gives the same differences (|Δe|/e = 2.0–2.2×10⁻³, 0.30–0.40% less energy above freeze-out, a 0.19% smaller freeze-out volume, 0.18–0.22% fewer charged hadrons; flow and ⟨pT⟩ unchanged), so the offset is a property of the GPU path, not of one backend; see the M3 Max section of docs/MUSIC_CPU_vs_GPU.md.

What this means: the GPU back-ends can be used for physics production. Keep a campaign on one path: the 0.2% offset in multiplicity is negligible within a campaign but does not average out when GPU and CPU samples are mixed, and name the path when comparing with CPU MUSIC results.

Still missing: a real large-scale campaign compared with the CPU build at the level of final-state observables (thousands of events, several centralities and systems; the test above is 6 events, 0–10% Au+Au, one jet pT̂ window). GPU campaigns of ~1000 events have run, but not against a CPU counterpart. Also open: a GPU without code for its architecture (e.g. a V100 with a build lacking sm_70) still runs on with meaningless fields until the CUDA fail-fast (PR #14) is merged; check the logs for "no kernel image".

Citation

Publishing results computed with MUSIC

New results computed with MUSIC may necessitate modifying certain parts of the code (the routines compute hadronic observables, transport coefficients, initial conditions, ...). To avoid any confusion about the physics behind MUSIC, we ask that publications containing results computed with MUSIC specify the initial conditions, equation of state, transport coefficients (first and second order), freeze-out criteria and other information necessary to reproduce the results.

The following papers should be cited when referring to MUSIC:

1. Schenke, S. Jeon, C. Gale. "3+1D hydrodynamic simulation of relativistic heavy-ion collisions" Phys.Rev.C 82, 014903 (2010) [arXiv:1004.1408]
2. Schenke, S. Jeon, C. Gale. "Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics” Phys.Rev.Lett. 106, 042301 (2011) [arXiv:1009.3244]
3. J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale. "Production of photons in relativistic heavy-ion collisions” Phys. Rev. C 93, 044906 (2016) [arXiv:1509.06738]

If simulations include net baryon diffusion, please cite in addition,

• G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, "Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions," Phys. Rev. C98, 034916 (2018) arXiv:1804.10557 [nucl-th]

The freeze-out surface at constant temperature or energy density is generated by Cornelius. When using this feature, please cite

• P. Huovinen and H. Petersen, "Particlization in hybrid models," Eur. Phys. J. A48, 171 (2012), arXiv:1206.3371 [nucl-th]

Modifying MUSIC

MUSIC users may need to modify various parts of the code. The GPL license allows for such modifications to be made and distributed (subject to restrictions that can be found in the license). Nevertheless, we ask and require that you state clearly any scientifically relevant modifications made to the code when you share this modified version of MUSIC. This is to avoid any misunderstanding about the physics behind the code.

Contributing to MUSIC

Should you make changes to MUSIC that you believe would benefit from being shared with the wider community, please get in touch with the code's maintainers through the MUSIC website. Contributions will be considered whenever possible. While no promises can be made regarding the inclusion of these modifications in the main version of the code, they could be made available separately on the MUSIC website for other users' convenience. Bug fixes will be considered separately and will be applied promptly.

Compilation

The MUSIC code can be compiled using standard cmake.

Alternatively, you can compile the MUSIC code using a makefile. In this way, please make sure the information about the openMP and the directory of the GSL library is correct in the src/GNUmakefile. Then one can compile the code by typing make.

Run MUSIC with an input file

An input file is required that contains the line EndOfData, preceded by a list of parameter names and values, one per line, with parameter names and values separated by a space or tab. If omitted, each parameter will be assigned a default value. Example input files can be found under the folder example_inputfiles/.

The list of possible parameters and default values, along with a brief description, can be found in utilities/music_parameters_dict.py.

To run hydro simulation, one can start by running the python script generate_music_inputfile.py (can be found in utilites folder) to generate input files and job running script. The job script can run under bash and qsub systems. For some help information about generate_music_inputfile.py, one can simply type:

generate_music_inputfile.py -h

Run MUSIC on multiple CPU cores

MUSIC uses OpenMP for parallelization. For example, to run on two processors and use the sample input file, type:

export OMP_NUM_THREADS=2
./MUSIChydro input_example

macOS (AppleClang) note. AppleClang ships without an OpenMP runtime. Install Homebrew's libomp keg once:

brew install libomp

CMake auto-detects it via brew --prefix libomp and links the resulting binary against /opt/homebrew/opt/libomp/lib/libomp.dylib. If libomp is absent the build falls back to serial execution and prints a notice during configure. Once the prerequisites are installed, you can build the package using:

make -j 10 #Adjust 10 to the number of cores available.

The result will be an executable named MUSIChydro.

---

GPU Acceleration

MUSIC 3.1's hydrodynamic evolution has been ported to run on the GPU, with two interchangeable back-ends selected at build time. No new physics is introduced — these are ports of the existing CPU algorithms (see the disclaimer at the top of this file).

Back-end	Target hardware	Build flag	Documentation
CUDA	NVIDIA GPUs (discrete + Grace-Blackwell)	-DUSE_CUDA=ON	README_CUDA.md
Metal	Apple Silicon (M1 or later)	-DUSE_METAL=ON	README_Metal.md

Both back-ends share the same seven compute kernels and the same host dispatch logic in src/advance.cpp (the CUDA back-end is a direct port of the Metal one). The plain CPU build (no flag) is unchanged and is the reference the GPU paths are validated against.

Kokkos port (separate branch, behind). A performance-portable Kokkos back-end (one kernel source for NVIDIA, AMD and Intel GPUs and multicore CPUs) exists on branch XSCAPE-KoKKos (PR #8; see README_KoKKos.md there). On the GB10 it reproduces the CPU eps_max evolution as closely as native CUDA (6.4e-4), at ~0.8× native CUDA's speed. HIP and SYCL are wired up but untested. It is not merged into XSCAPE and branched off before the X-SCAPE production work, so it lacks the bug fixes and features the CUDA and Metal back-ends have since received:

• the GPU fix for stalled evolutions in dilute regions (vacuum cells at rest, non-finite W/Π guard; PR #13, docs/VacReset_BUG.md);
• the parallel, deterministic freeze-out surface with the surface pressure initialised (#12);
• the X-SCAPE jet source slot (liquefier droplets), the freeze_out_surface switch (#9), the source-fill speed-ups (#10), the StringFind4 EOF guard (#11), releasing the evolution store's memory (#15), less host work per run (#19), and MusicCheckpoint (branch music2_restart);
• the CPU-vs-GPU validation above (#16).

Use CUDA or Metal for productions; the Kokkos branch needs these changes merged in, and its own validation repeated, before it can be used for one.

What is ported to the GPU

When the configuration falls inside the supported matrix, the entire per-cell hydro timestep runs on the GPU (the per-cell CPU loop is skipped):

• Ideal evolution — KT flux (MakeDeltaQI) + conserved→primitive Newton reconstruction
• Viscous source divergence (MakeWSource)
• Shear stress π^{μν} — full second-order Israel–Stewart update
• Bulk pressure Π
• Velocity gradients θ, aᵘ, σ^{μν}
• T-dependent η/s (modes 0, 1, 2, 3, 11) and ζ/s (modes 0, 1, 2, 3, 7, 8, 9, 10)
• Second-order coupling terms; QuestRevert regulator
• Hydro source terms (source evaluated on CPU, integrated on GPU)
• Zero net-baryon density (rhob = 0 EOS); 2D boost-invariant and full 3+1D grids

What stays on the CPU (not ported)

These configurations transparently fall back to the CPU path — correct, just slower:

• Baryon diffusion (turn_on_diff = 1)
• Vorticity terms (include_vorticity_terms = 1)
• Finite net-baryon EOS / μ_B-dependent shear (the GPU EOS table is sampled at rhob = 0)
• μ_B-dependent shear multiplier (muB_dependent_shear_to_s != 0)
• Cooper–Frye freeze-out / Cornelius surface finding, initial-condition construction, evolution output

Settings for X-SCAPE productions

The rest of a run's host work, and how to switch it, for the X-SCAPE productions. These changes came with the bit-identical CPU fixes of js-contrib's docs/Plans/REPORT_CPU_GPU.md, PR #19 and branch music2_restart. Defaults keep the earlier behaviour.

setting	default	effect
output_momentum_anisotropy_vs_etas (input file, set_parameter)	1	0 skips Cell_info::output_momentum_anisotropy_vs_etas at the four early steps of a run: the momentum_anisotropy_tau_*.dat / eccentricities_evo_*_tau_*.dat files, a serial pass over every cell and a GPU → host copy each. X-SCAPE: <Hydro><MUSIC><output_momentum_anisotropy_vs_etas>; the js-contrib production XMLs set 0.
late_delta_tau, late_tau_switch (input file, set_parameter)	0, 0	A larger time step from late_tau_switch on: late_delta_tau [fm/c] replaces Delta_Tau, and facTau and output_evolution_every_N_timesteps are scaled with it so that the freeze-out checks and the stored frames keep their spacing in fm/c (0.02 → 0.05 with 5 → 2 keeps 0.1 fm/c). The switch waits for a step that is a frame step before and after it; the first check after it gets its true spacing; the sources are told (HydroSourceBase::hydro_dtau_changed, X-SCAPE's jet source follows it); DATA is restored when the run ends. 0: a constant step. A checkpoint after the switch is refused. X-SCAPE: <Hydro><MUSIC><late_delta_tau>, <late_tau_switch>. Not bit-identical: the time discretization changes (README_CUDA.md, Phase 7).
store_hydro_info_frames (input file, set_parameter)	1	With store_hydro_info_in_memory 1, 0 keeps the store's grid and τ bookkeeping (get_ntau, hydro_tau_max) but packs and stores no frame and reserves no memory for one: for an instance whose evolution nobody reads. X-SCAPE: MpiMusic::set_store_evolution(false).
MUSIC::set_checkpoint_tau(tau)	0 (off)	Every run keeps a MusicCheckpoint (src/music_checkpoint.h) at its last freeze-out step whose RK step ends by τ, after that step's freeze-out stage: snap_prev/snap_curr (a copy on the GPU), the surface so far, the stop-test bookkeeping, the store's frames, and a key (the parameters that shape the evolution, a hash of the initial arena, the initial-state source's fingerprint). get_checkpoint() returns it.
MUSIC::set_resume_from(checkpoint)	—	The next run starts from the checkpoint and enters its step at AdvanceRK, if the key matches, the needed surface and frames are there, and its own jet source adds nothing before that step. Otherwise it runs from τ0. get_resume_status() says which. One run only. X-SCAPE: MpiMusic::set_checkpoint_tau / set_restart_from, the jet leg of the two-stage chain starting from the background's state. GPU runs only.
MUSIC_RESTART_VERIFY=1 (env)	unset	A run given a checkpoint runs from τ0 instead, and at the checkpoint's step compares its state with it bit for bit (GPU state, surface, max e and counters, store frames); the log says identical to the checkpoint or what differs.
MUSIC_FREEZEOUT_GATHER=0 (env)	unset	The freeze-out search copies every field of every cell of the checked step from the GPU, instead of e alone and the fields of the cubes the surface crosses (below). The surface is the same either way; for comparisons.

Changes that need no setting:

• Charge lists that are not evolved (turn_on_rhob, turn_on_QS) no longer trigger a string source fill of zeros.
• The stored frames are packed straight into the end of the in-memory store, whose reserved memory is madvise(MADV_HUGEPAGE)'d on Linux.
• On discrete CUDA GPUs the pack goes through a device buffer and a pinned host buffer instead of managed memory. MUSIC_CUDA_FORCE_DISCRETE=1 exercises that path on a coherent GPU.
• The no-surface stop test's edge check is a GPU reduction (CUDA kernel; a host scan of the shared buffer on Metal), with no grid copy.
• The freeze-out search (Cornelius, every facTau steps) no longer copies the GPU state into the host arenas and then both arenas into the "previous step" buffers. Two buffers take the checked steps in turn (FreezeoutSteps), and while the GPU holds the state only e comes to the host: the cubes the surface crosses are found from e, and the other fields are gathered for their corner cells only, from the GPU state and from a copy of the previously checked step kept on the GPU (CUDA kernel gpu_gather_cells; a host gather from the shared buffers on Metal). The surface is bit-identical. Not for the boost-invariant search or output_vorticity, which read more than the cube corners: they take full copies, into the same two buffers.

The Metal parts of the checkpoint and of the edge check were built and checked on an M3 Max (2026-10-11): outputs bit-identical to the code before these changes.

Benchmarks (indicative)

With (more relevant for discete GPUs):

output_evolution_every_N_timesteps>1</output_evolution_every_N_timesteps>

Back-end	GPU	CPU baseline	Speedup @ single O+O event tauMax > 3fm
Metal	Apple M3 Max	16-thread	~8× (wall)
CUDA	NVIDIA GB10	20-thread	~5.5× (wall)
	RTX3090	48-thread	~3.5× (wall)

These are not a controlled head-to-head: the two back-ends run on different machines, different CPUs/thread counts, so read each number as a speedup over its own host CPU, not as Metal-vs-CUDA. The "wall" figures are total run time (including ~0.5 s of one-time GPU init); Full per-phase and per-grid tables, the methodology, and the reproduction scripts (tests/cuda_vs_cpu_bench*.sh, tests/cuda_perstep_bench.sh) are in the per-back-end READMEs.

See the per-back-end READMEs (CUDA, Metal) for build instructions, kernel-level details, benchmarks, and the full support matrix with what each missing piece would take to port.