PLUTO Hydrodynamical Code

• PLUTO hydrodynamical code is a modular, grid-based solver that employs high-resolution shock-capturing finite-volume methods for HD/MHD simulations in various coordinate systems.
• It integrates specialized modules for radiation, self-gravity, and chemistry, enabling applications from accretion discs and protoplanetary systems to radiative shocks and photoionized flows.
• Benchmarking studies and extensions like gPLUTO and PyPLUTO underscore its robust architecture, flexible multi-physics capabilities, and adaptability to modern computational techniques.

PLUTO is a modular, grid-based code for computational astrophysics that solves the equations of hydrodynamics and magnetohydrodynamics using shock-capturing Godunov-type schemes on structured grids in Cartesian, cylindrical, and spherical coordinates. In the literature represented here, it appears both as a general-purpose finite-volume engine for compressible HD/MHD and as a platform to which radiation hydrodynamics, photoionization, self-gravity, and post-processing tools have been added for applications ranging from thin accretion discs and protoplanetary disks to radiative shocks, H II regions, and self-gravitating collapse (Kolb et al., 2013, Parthasarathy et al., 2014, Mattia et al., 16 Jan 2025).

1. Core formulation

In the hydrodynamical configuration used for thin-disc simulations, PLUTO is described as a grid-based, Godunov-type, high-resolution shock-capturing finite-volume code for solving the hydrodynamic equations in the presence of gravity (Parthasarathy et al., 2014). Its generic conservation-law form is written as

$$\frac{\partial U}{\partial t} = - \nabla \cdot T(U) + S(U),$$

where $U$ is the vector of conserved variables, $T(U)$ is the flux tensor, and $S(U)$ contains source terms such as gravity. In the HD module, the governing equations are the Euler equations with gravity in conservative form: mass conservation,

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0,$$

momentum conservation,

$$\frac{\partial (\rho \mathbf{v})}{\partial t} + \nabla \cdot (\rho \mathbf{v}\mathbf{v} + p \mathbf{I}) = -\rho \nabla \Phi,$$

and total-energy conservation with gravitational work (Parthasarathy et al., 2014).

A complementary description in the radiation-module literature presents PLUTO as a 3D, finite-volume, grid-based code that supports hydrodynamics and magnetohydrodynamics, with source terms for gravity and other physics written in the same conservative framework (Kolb et al., 2013). Across the papers considered here, closure depends on the physical problem. Ideal-gas equations of state are used in hydrodynamical and radiation-hydrodynamical studies, while an isothermal equation of state,

$P = C_s^2 \rho,$

is used in comparative HD/MHD tests designed to align PLUTO with pseudo-spectral reference solvers (Biswas et al., 2024).

This formulation makes PLUTO a finite-volume solver centered on conservation, Riemann problems at cell interfaces, and explicit treatment of hyperbolic transport. A plausible implication is that the code’s identity is defined less by a single physical regime than by a reusable conservative infrastructure into which different source terms, closures, and physics modules can be inserted.

2. Numerical architecture and discretization

PLUTO’s numerical design is repeatedly characterized as multi-physics, multi-algorithm, and high-resolution. The code supports Cartesian, cylindrical, and spherical coordinates, and the radiation-module literature further notes support for non-equidistant grids (Kolb et al., 2013, Biswas et al., 2024). Its finite-volume update uses reconstructed interface states, approximate Riemann solvers, and explicit time integration; in MHD configurations it also supports constrained transport for the solenoidal constraint on $\mathbf{B}$ (Biswas et al., 2024).

Different studies emphasize different solver choices, reflecting PLUTO’s configurability rather than a single canonical setup.

Study	Physics and geometry	Example numerical choices
Thin accretion discs (Parthasarathy et al., 2014)	HD, 2.5D spherical $(R,\theta)$	127 uniformly spaced points in each direction
GMHD3D comparison (Biswas et al., 2024)	Isothermal HD/MHD, Cartesian, periodic	PARABOLIC reconstruction, HLLD, RK3, constrained transport, explicit viscosity and resistivity
AGN jet comparison (Maragkakis et al., 21 Oct 2025)	3D Cartesian HD	HLLC, linear reconstruction, RK2; shock flattening with HLL and MINMOD in strong shocks

The thin-disc study emphasizes the HRSC framework itself rather than a specific solver stack: PLUTO computes interface fluxes via approximate Riemann solvers, handles shocks and discontinuities while smearing them over only a few grid cells, and uses a structured grid in spherical coordinates (Parthasarathy et al., 2014). The GMHD3D comparison provides a more explicit example of a high-order PLUTO4.4 configuration: finite volume with PARABOLIC reconstruction, HLLD for isothermal MHD, RK3 time stepping, and constrained transport with several CT_EMF_AVG choices that materially affect magnetic-field evolution and dissipation (Biswas et al., 2024). The AGN-jet comparison, by contrast, uses PLUTO in a purely hydrodynamic, non-relativistic, static-grid mode with HLLC, second-order linear reconstruction, RK2, and a localized switch to the more diffusive HLL solver plus a MINMOD limiter in strong shocks (Maragkakis et al., 21 Oct 2025).

These examples show that PLUTO’s “hydrodynamical code” identity is inseparable from configurability. The same code base can be deployed as a shock-robust second-order Eulerian solver, as a higher-order isothermal MHD code with constrained transport, or as the transport backbone for much richer radiation and self-gravity systems.

3. Radiation-hydrodynamics extensions

A major strand of PLUTO development concerns radiation hydrodynamics. One line of work adds a grey, two-temperature flux-limited diffusion solver in the co-moving frame, with hydrodynamics handled by PLUTO’s Godunov schemes and radiation transport handled by a fully implicit scheme whose linear systems are solved either by SOR or by PETSc-based matrix solvers (Kolb et al., 2013). In that implementation, the module supports Cartesian, cylindrical, and spherical coordinates, non-equidistant grids, and parallel execution with MPI, and it was validated on linear-diffusion tests, gas–radiation coupling tests, radiative shocks, and accretion-disc convection (Kolb et al., 2013).

A second line of work introduces a gray, two-moment radiation-hydrodynamics module based on the M1 closure, designed as a PLUTO module for simulations of planet formation in circumstellar disks (Fuksman et al., 2020). The radiation subsystem evolves $E_r$ and $U$0 and uses operator splitting, substepping, and the reduced speed of light approximation. Radiation transport is advanced explicitly with Godunov-type solvers, while radiation–matter interaction terms are integrated with IMEX schemes. The module is formulated for 1D, 2D, and 3D in Cartesian, spherical, and cylindrical coordinates, and can be implemented on rotating frames (Fuksman et al., 2020). This is technically distinct from the FLD approach: FLD reduces transport to a diffusion problem with a flux limiter, whereas the M1 closure retains a two-moment hyperbolic radiation system and is explicitly designed to cover both the diffusion and free-streaming limits.

PLUTO has also been extended from two-temperature gas–radiation coupling to a three-temperature scheme that evolves gas, dust, and radiation as distinct thermal components (Muley et al., 2023). In that work, the code is the numerical backbone for a “three-temperature” radiation-hydrodynamics module tailored to protoplanetary disks around T Tauri stars. The formulation includes separate energy equations for gas and dust, a grey radiation field with M1 closure, gas–dust collisional coupling, operator splitting, IMEX radiation stepping, and a local Newton–Raphson solve for stiff energy exchange terms (Muley et al., 2023). The resulting framework was exercised in 0D, 1D, 2D, and 3D tests, including dusty radiative shocks, disk self-shadowing, and planet-driven spirals.

In the relativistic domain, PLUTO also hosts a grey relativistic radiation hydrodynamics and magnetohydrodynamics module with M1 closure (Fuksman et al., 2019). There, the radiation transfer equations are solved under the grey approximation and integrated with an implicit-explicit scheme in which radiation–matter interaction terms are implicit and transport is handled by the same Godunov-type solvers used by PLUTO’s HD/RMHD infrastructure. The paper introduces a radiation HLLC solver for optically thin transport and demonstrates applicability in Cartesian, spherical, and cylindrical geometries, on single processors, parallel architectures, and adaptive grids via CHOMBO (Fuksman et al., 2019).

Taken together, these modules show that PLUTO’s hydrodynamical core has repeatedly served as the transport layer for increasingly elaborate radiation closures rather than being replaced by them. This suggests that the finite-volume HD/MHD engine is the stable architectural element, while radiation physics is layered on top through operator splitting and implicit local or global solves.

4. Photoionization, chemistry, and coupled microphysics

Another branch of development couples PLUTO to non-equilibrium ionization, photoionization, and external plasma microphysics. A general-purpose time-dependent ionisation network and a radiation transport module were introduced in the MHD code PLUTO for temperatures ranging from $U$1 to $U$2 K, including all ionisation states of H, He, C, N, O, Ne, Mg, Si, S, and Fe (Sarkar et al., 2020). In that implementation, radiation loss for each ion-electron pair is calculated using CLOUDY-17 data on-the-fly, photo-ionisation and charge exchange are the chemical heating mechanisms, and the radiation transport module employs short characteristics assuming spherical symmetry (Sarkar et al., 2020). The ionisation network is compatible with full 3D, whereas the radiation module requires spherical symmetry.

A separate route to detailed microphysics is the PLUTO–CLOUDY Interface, or TPCI (Salz et al., 2015). TPCI couples PLUTO’s hydrodynamics or MHD to CLOUDY’s equilibrium photoionization and spectral-synthesis machinery, constructing a photoionization hydrodynamics solver aimed at strongly irradiated photoevaporative flows. In the iterative workflow described there, CLOUDY first solves the equilibrium state of the medium for a given radiation field and returns a net radiative heating or cooling rate; in the second step, that source term influences the hydrodynamic evolution calculated by PLUTO (Salz et al., 2015). The one-dimensional implementation was validated on photoevaporation of a cool hydrogen cloud, cooling of coronal plasma, formation of a Stroemgren sphere, and the evaporating atmosphere of a hot Jupiter (Salz et al., 2015).

The MAKEMAKE and SEDNA framework extends PLUTO in yet another direction by combining continuum radiation transport and hydrogen photoionization on static one-dimensional, two-dimensional, and three-dimensional grids in Cartesian, cylindrical, and spherical coordinates (Kuiper et al., 2020). Each module splits the radiation field into a direct component from a point source, solved with ray tracing, and a diffuse component, solved with a three-dimensional FLD solver. For continuum radiation, the FLD solver supports either an equilibrium one-temperature approach or a linearization two-temperature approach; for photoionization, the diffuse EUV field from direct recombination into hydrogen’s ground state is evolved as an alternative to the on-the-spot approximation (Kuiper et al., 2020). Radiation and ionization forces are then passed into PLUTO through the external-acceleration interface.

These systems demonstrate that PLUTO has functioned not only as a standalone hydrodynamical code but also as a coupling hub for external or semi-external microphysics solvers. This suggests a division of labor in which PLUTO supplies conservative fluid transport, geometry handling, and parallel domain decomposition, while specialized modules supply thermochemistry, opacity physics, or photon transport.

5. Self-gravity and elliptic solvers

Self-gravity has become a major part of the PLUTO ecosystem. A first dedicated self-gravity module implements an iterative solver for Poisson’s equation based on parabolic relaxation and Runge–Kutta–Legendre super-time-stepping, embedded in a multigrid framework (Mandal et al., 2023). The solver is described as fully parallel and rapidly convergent, and the authors present it as a practical alternative to conventional iterative solvers such as Gauss–Seidel and successive over-relaxation (Mandal et al., 2023). On structured grids, the scheme uses second-order cell-centered discretization of the Laplacian, restriction and prolongation between multigrid levels, and a fourth-order centered finite difference for gravitational acceleration derived from the potential. It is coupled to PLUTO’s Euler update either once per timestep or at each RK stage (Mandal et al., 2023).

That self-gravity module is validated on a broad range of standard self-gravitating astrophysical problems, including a 3D analytical model potential, Jeans instability, Evrard collapse, non-rotating isothermal collapse, rotating collapse, and the Boss–Bodenheimer–Burkert fragmentation test (Mandal et al., 2023). The reported behavior includes second-order convergence for the potential, excellent agreement with analytic oscillation and collapse times in the Jeans problem, total-energy conservation to $U$3 in Evrard collapse, and successful formation of a centrifugally supported disk without an artificial ring in the rotating-cloud test (Mandal et al., 2023).

A later development targets non-Cartesian self-gravity directly. A fast spectral-multigrid Poisson solver was implemented in PLUTO for spherical and cylindrical geometries with either uniform or logarithmic radial discretizations (Mandal et al., 16 Jun 2026). The method uses azimuthal Fourier decomposition to turn the 3D Poisson equation into a set of 2D Helmholtz equations, which are then solved with a geometrically consistent multigrid algorithm. Vacuum boundary conditions are implemented through a screening-mass approach, and the solver is validated against analytical solutions and dynamical test problems in both spherical and cylindrical geometries (Mandal et al., 16 Jun 2026). The paper reports second-order convergence, excellent agreement with reference solutions, and weak-scaling tests up to 4096 cores in which the Poisson solve remains subdominant to the magnetohydrodynamic update cost (Mandal et al., 16 Jun 2026).

The coexistence of these two gravity strategies is significant. The earlier module establishes a general self-gravity capability on structured grids; the later one specializes to the non-Cartesian regimes that are particularly important for disks, collapse, and strongly non-axisymmetric structures. A plausible implication is that PLUTO’s gravity infrastructure is evolving from “self-gravity added to a fluid code” toward geometry-aware elliptic solvers designed around the same curvilinear grids used by the hydrodynamics.

6. Benchmarks, analysis ecosystem, and computational trajectory

PLUTO has been evaluated in multiple code-comparison settings. In the exoALMA benchmark of forward modeling for protoplanetary disks, PLUTO was one of four grid-based hydrodynamics codes used to simulate a disk with an embedded giant planet (Bae et al., 25 Apr 2025). That study found strong consistency between the hydrodynamics codes, particularly in the density and velocity perturbations associated with planet-driven spirals, and also found good agreement between the two radiative-transfer codes used in post-processing: the disk temperature in mcfost and RADMC-3D agreed within $U$4 everywhere in the domain, brightness-temperature differences in synthetic $U$5CO channel maps stayed within $U$6 K, retrieved planet radial and azimuthal locations showed only a few percent of scatter, and velocity perturbations varied $U$7 among the models (Bae et al., 25 Apr 2025). The abstract concludes that any combination of the tested hydrodynamics and radiative-transfer codes can be used to reliably model and interpret planet-driven kinematic perturbations (Bae et al., 25 Apr 2025).

A different comparison, between PLUTO4.4 and the pseudo-spectral multi-GPU solver GMHD3D, reaches a more differentiated conclusion (Biswas et al., 2024). There, PLUTO4.4 and GMHD3D are reported to produce remarkably similar results for hydrodynamic and magnetohydrodynamic problems such as instability growth, three-dimensional turbulent dynamics, oscillation of kinetic and magnetic energy, and recurrence dynamics. At the same time, the authors state that the pseudo spectral solver GMHD3D is significantly more superior than the grid based solver PLUTO4.4 for certain category of physics problems, especially in periodic-box settings where higher spatial order and lower numerical diffusion dominate the comparison (Biswas et al., 2024). This is an important corrective to the common misconception that a general-purpose astrophysical finite-volume code should always be judged by the same criteria as a spectral turbulence code.

Post-processing has its own dedicated companion in the form of PyPLUTO, a Python package tailored for efficient loading, manipulation, and visualization of outputs produced with the PLUTO code (Mattia et al., 16 Jan 2025). PyPLUTO uses memory mapping to optimize data loading, supports binary, HDF5, VTK, and ASCII outputs, includes dedicated classes for fluid data and particle data, provides FFT and differential-operator utilities, and offers both scripted visualization and a PyQt6-based graphical user interface for interactive 1D and 2D exploration (Mattia et al., 16 Jan 2025). Its scope includes the particle modules of PLUTO, allowing cosmic rays, Lagrangian particles, and dust particles to be loaded and analyzed alongside the fluid fields (Mattia et al., 16 Jan 2025).

At the implementation level, the code base is also undergoing architectural change. The GPU-oriented successor gPLUTO is presented as a complete rewrite in C++ that uses OpenACC to accelerate PLUTO on NVIDIA GPUs while retaining the finite-volume formulation for Eulerian MHD methods in multiple spatial dimensions (Rossazza et al., 25 Nov 2025). The paper frames these results as a first look and emphasizes that the GPU implementation is not a minor port but a new internal realization of the PLUTO approach (Rossazza et al., 25 Nov 2025).

The cumulative picture is of a code that has been repeatedly benchmarked, extended, analyzed, and reimplemented without abandoning its underlying finite-volume, Godunov-type identity. That continuity is central to understanding PLUTO: it is not a single solver frozen at one stage of astrophysical numerics, but a persistent hydrodynamical framework around which new radiation, gravity, chemistry, analysis, and hardware layers continue to be built.