
We describe an implementation of a particle physics module available for the
PLUTO code, appropriate for the dynamical evolution of a plasma consisting of a
thermal fluid and a nonthermal component represented by relativistic charged
particles, or cosmic rays (CR). While the fluid is approached using standard
numerical schemes for magnetohydrodynamics, CR particles are treated
kinetically using conventional ParticleInCell (PIC) techniques.
The module can be used to describe either test particles motion in the fluid
electromagnetic field or to solve the fully coupled MHDPIC system of equations
with particle backreaction on the fluid as originally introduced by
\cite{Bai_etal.2015}. Particle backreaction on the fluid is included in the
form of momentumenergy feedback and by introducing the CRinduced Hall term in
Ohm's law. The hybrid MHDPIC module can be employed to study CR kinetic
effects on scales larger than the (ion) skin depth provided the Larmor gyration
scale is properly resolved. When applicable, this formulation avoids to resolve
microscopic scales offering a substantial computational saving with respect to
PIC simulations.
We present a fullyconservative formulation which is secondorder accurate in
time and space and extends to either RungeKutta (RK) or
cornertransportupwind (CTU) timestepping schemes (for the fluid) while a
standard Boris integrator is employed for the particles. For highlyenergetic
relativistic CRs and in order to overcome the time step restriction a novel
subcycling strategy that retains secondorder accuracy in time is presented.
Numerical benchmarks and applications including Bell instability, diffusive
shock acceleration and test particle acceleration in reconnecting layers are
discussed.

We present the results of stateoftheart simulations of recollimation
shocks induced by the interaction of a relativistic jet with an external
medium, including the effect of radiative losses of the shocked gas. Our
simulations confirm that  as suggested by earlier semianalytical models 
the postshock pressure loss induced by radiative losses may lead to a
stationary equilibrium state characterized by a very strong focusing of the
flow, with the formation of quite narrow nozzles, with crosssectional radii as
small as $10^{3}$ times the length scale of the jet. We also study the
timedependent evolution of the jet structure induced of a density perturbation
injected at the flow base. The setup and the results of the simulations are
particularly relevant for the interpretation of the observed rapid variability
of the $\gamma$ray emission associated to flat spectrum radio quasars. In
particular, the combined effects of jet focusing and Doppler beaming of the
observed radiation make it possible to explain the subhour flaring events such
as that observed in the FSRQ PKS 1222+216 by MAGIC.

We investigate magnetohydrodynamic turbulence driven by the magnetorotational
instability (MRI) in Keplerian disks with a nonzero net azimuthal magnetic
field using shearing box simulations. As distinct from most previous studies,
we analyze turbulence dynamics in Fourier (${\bf k}$) space to understand its
sustenance. The linear growth of MRI with azimuthal field has a transient
character and is anisotropic in Fourier space, leading to anisotropy of
nonlinear processes in Fourier space. As a result, the main nonlinear process
appears to be a new type of angular redistribution of modes in Fourier space 
the \emph{nonlinear transverse cascade}  rather than usual direct/inverse
cascade. We demonstrate that the turbulence is sustained by interplay of the
linear transient growth of MRI (which is the only energy supply for the
turbulence) and the transverse cascade. These two processes operate at large
length scales, comparable to box size and the corresponding small wavenumber
area, called \emph{vital area} in Fourier space is crucial for the sustenance,
while outside the vital area direct cascade dominates. The interplay of the
linear and nonlinear processes in Fourier space is generally too intertwined
for a vivid schematization. Nevertheless, we reveal the \emph{basic subcycle}
of the sustenance that clearly shows synergy of these processes in the
selforganization of the magnetized flow system. This synergy is quite robust
and persists for the considered different aspect ratios of the simulation
boxes. The spectral characteristics of the dynamical processes in these boxes
are qualitatively similar, indicating the universality of the sustenance
mechanism of the MRIturbulence.

A significant fraction of extended radio sources presents a peculiar Xshaped
radio morphology: in addition to the classical double lobed structure, radio
emission is also observed along a second axis of simmetry in the form of
diffuse wings or tails. In a previous investigation we showed the existence of
a connection between the radio morphology and the properties of the host
galaxies. Motivated by this connection we performed twodimensional numerical
simulations showing that Xshaped radio sources may naturally form as a jet
propagates along the major axis a highly elliptical density distribution,
because of the fast expansion of the cocoon along the minor axis of the
distribution.
We intend to extend our analysis by performing threedimensional numerical
simulations and investigating the role of different parameters is determining
the formation of the Xshaped morphology.
The problem is addressed by numerical means, carrying out threedimensional
relativistic magnetohydrodynamic simulations of bidirectional jets propagating
in a triaxial density distribution.
We show that only jets with power $\lesssim 10^{44}$ erg s$^{1}$ can give
origin to an Xshaped morphology and that a misalignment of $30^o$ between the
jet axis and the major axis of the density distribution is still favourable to
the formation of this kind of morphology. In addition we compute synthetic
radio emission maps and polarization maps.
In our scenario for the formation of Xshaped radio sources only low power
FRII can give origin to such kind of morphology. Our synthetic emission maps
show that the different observed morphologies of Xshaped sources can be the
result of similar structures viewed under different perspectives.

We examine the relationship between magnetic flux generation, taken as an
indicator of largescale dynamo action, and magnetic helicity, computed as an
integral over the dynamo volume, in a simple dynamo. We consider dynamo action
driven by MagnetoRotational Turbulence (MRT) within the shearingbox
approximation. We consider magnetically open boundary conditions that allow a
flux of helicity in or out of the computational domain. We circumvent the
problem of the lack of gauge invariance in open domains by choosing a
particular gauge  the winding gauge  that provides a natural interpretation
in terms of average winding number of pairwise field lines. We use this gauge
precisely to define and measure the helicity and helicity flux for several
realizations of dynamo action. We find in these cases, that the system as a
whole does not break reflectional symmetry and the total helicity remains small
even in cases when substantial magnetic flux is generated. We find no
particular connection between the generation of magnetic flux and the helicity
or the helicity flux through the boundaries. We suggest that this result may be
due to the essentially nonlinear nature of the dynamo processes in MRT.

Extragalactic radio sources have been classified into two classes,
FanaroffRiley I and II, which differ in morphology and radio power. Strongly
emitting sources belong to the edgebrightened FR II class, and weakly emitting
sources to the edgedarkened FR I class. The origin of this dichotomy is not
yet fully understood. Numerical simulations are successful in generating FR II
morphologies, but they fail to reproduce the diffuse structure of FR Is.
By means of hydrodynamical 3D simulations of supersonic jets, we investigate
how the displayed morphologies depend on the jet parameters. Bow shocks and
Mach disks at the jet head, which are probably responsible for the hot spots in
the FR II sources, disappear for a jet kinetic power L_kin < 10^43 erg/s. This
threshold compares favorably with the luminosity at which the FR I/FR II
transition is observed.
The problem is addressed by numerical means carrying out 3D HD simulations of
supersonic jets that propagate in a nonhomogeneous medium with the ambient
temperature that increases with distance from the jet origin, which maintains
constant pressure.
The jet energy in the lower power sources, instead of being deposited at the
terminal shock, is gradually dissipated by the turbulence. The jets spread out
while propagating, and they smoothly decelerate while mixing with the ambient
medium and produce the plumes characteristic of FR I objects.
Threedimensionality is an essential ingredient to explore the FR I evolution
because the properties of turbulence in two and three dimensions are very
different, since there is no energy cascade to small scales in two dimensions,
and twodimensional simulations with the same parameters lead to FRIIlike
behavior.

We perform a linear stability analysis of magnetized rotating cylindrical jet
flows in the approximation of zero thermal pressure. We focus our analysis on
the effect of rotation on the current driven mode and on the unstable modes
introduced by rotation. We find that rotation has a stabilizing effect on the
current driven mode only for rotation velocities of the order of the Alfv\'en
velocity. Rotation introduces also a new unstable centrifugal buoyancy mode and
the "cold" magnetorotational instability. The first mode is analogous to the
Parker instability with the centrifugal force playing the role of effective
gravity. The magnetorotational instability can be present, but only in a very
limited region of the parameter space and is never dominant. The current driven
mode is characterized by large wavelenghts and is dominant at small values of
the rotational velocity, while the buoyancy mode becomes dominant as rotation
is increased and is characterized by small wavelenghts.

We present an approach to deriving global properties of accretion disks from
the knowledge of local solutions derived from numerical simulations based on
the shearing box approximation. The approach consists of a twostep procedure.
First a local solution valid for all values of the disk height is constructed
by piecing together an interior solution obtained numerically with an
analytical exterior radiative solution. The matching is obtained by assuming
hydrostatic balance and radiative equilibrium. Although in principle the
procedure can be carried out in general, it simplifies considerably when the
interior solution is fully convective. In these cases, the construction is
analogous to the derivation of the Hayashi tracks for protostars. The second
step consists of piecing together the local solutions at different radii to
obtain a global solution. Here we use the symmetry of the solutions with
respect to the defining dimensionless numbersin a way similar to the use of
homology relations in stellar structure theoryto obtain the scaling
properties of the various disk quantities with radius.

An Equation of State (\textit{EoS}) closes the set of fluid equations.
Although an ideal EoS with a constant \textit{adiabatic index} $\Gamma$ is the
preferred choice due to its simplistic implementation, many astrophysical fluid
simulations may benefit from a more sophisticated treatment that can account
for diverse chemical processes. Here, we first review the basic thermodynamic
principles of a gas mixture in terms of its thermal and caloric EoS by
including effects like ionization, dissociation as well as temperature
dependent degrees of freedom such as molecular vibrations and rotations. The
formulation is revisited in the context of plasmas that are either in
equilibrium conditions (local thermodynamic or collisional excitation
equilibria) or described by nonequilibrium chemistry coupled to optically thin
radiative cooling. We then present a numerical implementation of thermally
ideal gases obeying a more general caloric EoS with nonconstant adiabatic
index in Godunovtype numerical schemes.We discuss the necessary modifications
to the Riemann solver and to the conversion between total energy and pressure
(or viceversa) routinely invoked in Godunovtype schemes. We then present two
different approaches for computing the EoS.The first one employs rootfinder
methods and it is best suited for EoS in analytical form. The second one leans
on lookup table and interpolation and results in a more computationally
efficient approach although care must be taken to ensure thermodynamic
consistency. A number of selected benchmarks demonstrate that the employment of
a nonideal EoS can lead to important differences in the solution when the
temperature range is $50010^4$ K where dissociation and ionization occur. The
implementation of selected EoS introduces additional computational costs
although using lookup table methods can significantly reduce the overhead by a
factor $3\sim 4$.

We investigate the linear and nonlinear evolution of currentcarrying jets in
a periodic configuration by means of high resolution threedimensional
numerical simulations. The jets under consideration are strongly magnetized
with a variable pitch profile and initially in equilibrium under the action of
a forcefree magnetic field. The growth of currentdriven (CDI) and
KelvinHelmholtz (KHI) instabilities is quantified using three selected cases
corresponding to static, Alfvenic and superAlfvenic jets.
During the early stages, we observe largescale helical deformations of the
jet corresponding to the growth of the initially excited CDI mode. A direct
comparison between our simulation results and the analytical growth rates
obtained from linear theory reveals good agreement on condition that
highresolution and accurate discretization algorithms are employed.
After the initial linear phase, the jet structure is significantly altered
and, while slowlymoving jets show increasing helical deformations, larger
velocity shear are violently disrupted on a few Alfven crossing time leaving a
turbulent flow structure. Overall, kinetic and magnetic energies are quickly
dissipated into heat and during the saturated regime the jet momentum is
redistributed on a larger surface area with most of the jet mass travelling at
smaller velocities. The effectiveness of this process is regulated by the onset
of KHI instabilities taking place at the jet/ambient interface and can be held
responsible for vigorous jet braking and entrainment.

We consider the problem of convergence in stratified isothermal shearing
boxes with zero net magnetic flux. We present results with the highest
resolution todateup to 200 gridpoint per pressure scale heightthat show
no clear evidence of convergence. Rather, the Maxwell stresses continue to
decrease with increasing resolution. We propose some possible scenarios to
explain the lack of convergence based on multilayer dynamo systems.

We investigate linear dynamics of nonaxisymmetric perturbations in
incompressible, vertically stratified Keplerian discs with a weak vertical
magnetic field in the shearing box approximation. Perturbations are decomposed
into shearing waves whose evolution is followed via numerical integration of
the linearized ideal MHD equations. There are two basic modes in the system 
inertiagravity waves and magnetic mode, which displays the magnetorotational
instability (MRI). As distinct from previous studies, we introduce
`eigenvariables' characterizing each (counterpropagating) component of the
inertiagravity and magnetic modes, which are governed by a set of four first
order coupled ordinary differential equations. This allowed us to identify a
new process of linear coupling of the two above nonaxisymmetric modes due to
the disc's differential rotation. We did a comparative analysis of the dynamics
of nonaxisymmetric and axisymmetric magnetic mode perturbations. It is shown
that the growth of optimal and closetooptimal nonaxisymmetric harmonics of
this mode, having transient nature, can prevail over the exponential growth of
axisymmetric ones (i.e., over the axisymmetric MRI) during dynamical time. A
possible implication of this result for axisymmetric channel solutions is
discussed. Specifically, the formation of the channel may be affected/impeded
by nonaxisymmetric modes already at the early linear stage leading to its
untimely disruption  the outcome strongly depends on the amplitude and
spectrum of initial perturbation. So, this competition may result in an
uncertainty in the magnetic mode's nonlinear dynamics. It is also shown that a
maximum nonaxisymmetric growth is at vertical wavelengths close to the
scaleheight for which compressibility effects are important. This indirectly
suggests that compressibility plays a role in the dynamics of the
nonaxisymmetric MRI.

We perform a linear analysis of the stability of a magnetized relativistic
nonrotating cylindrical flow in the aproximation of zero thermal pressure,
considering only the m = 1 mode. We find that there are two modes of
instability: KelvinHelmholtz and current driven. The KelvinHelmholtz mode is
found at low magnetizations and its growth rate depends very weakly on the
pitch parameter. The current driven modes are found at high magnetizations and
the value of the growth rate and the wavenumber of the maximum increase as we
decrease the pitch parameter. In the relativistic regime the current driven
mode is splitted in two branches, the branch at high wavenumbers is
characterized by the eigenfunction concentrated in the jet core, the branch at
low wavenumbers is instead characterized by the eigenfunction that extends
outside the jet velocity shear region.

We present a numerical study of turbulence and dynamo action in stratified
shearing boxes with zero magnetic flux. We assume that the fluid obeys the
perfect gas law and has finite (constant) thermal diffusivity. We choose
radiative boundary conditions at the vertical boundaries in which the heat flux
is propor tional to the fourth power of the temperature. We compare the
results with the corresponding cases in which fixed temperature boundary
conditions are applied. The most notable result is that the formation of a
fully convective state in which the density is nearly constant as a function of
height and the heat is transported to the upper and lower boundaries by
overturning motions is robust and persists even in cases with radiative
boundary conditions. Interestingly, in the convective regime, although the
diffusive transport is negligible the mean stratification does not relax to an
adiabatic state.

We present a numerical study of turbulence and dynamo action in stratified
shearing boxes with zero mean magnetic flux. We assume that the fluid obeys the
perfect gas law and has finite (constant) thermal diffusivity. The calculations
begin from an isothermal state spanning three scale heights above and below the
midplane. After a long transient the layers settle to a stationary state in
which thermal losses out of the boundaries are balanced by dissipative heating.
We identify two regimes. A conductive regime in which the heat is transported
mostly by conduction and the density decreases with height. In the limit of
large thermal diffusivity this regime resembles the more familiar isothermal
case. Another, the convective regime, observed at smaller values of the thermal
diffusivity, in which the layer becomes unstable to overturning motions, the
heat is carried mostly by advection and the density becomes nearly constant
throughout the layer. In this latter constantdensity regime we observe
evidence for largescale dynamo action leading to a substantial increase in
transport efficiency relative to the conductive cases.

We present a description of the adaptive mesh refinement (AMR) implementation
of the PLUTO code for solving the equations of classical and special
relativistic magnetohydrodynamics (MHD and RMHD). The current release exploits,
in addition to the static grid version of the code, the distributed
infrastructure of the CHOMBO library for multidimensional parallel computations
over blockstructured, adaptively refined grids. We employ a conservative
finitevolume approach where primary flow quantities are discretized at the
cellcenter in a dimensionally unsplit fashion using the Corner Transport
Upwind (CTU) method. Time stepping relies on a characteristic tracing step
where piecewise parabolic method (PPM), weighted essentially nonoscillatory
(WENO) or slopelimited linear interpolation schemes can be handily adopted. A
characteristic decompositionfree version of the scheme is also illustrated.
The solenoidal condition of the magnetic field is enforced by augmenting the
equations with a generalized Lagrange multiplier (GLM) providing propagation
and damping of divergence errors through a mixed hyperbolic/parabolic explicit
cleaning step. Among the novel features, we describe an extension of the scheme
to include nonideal dissipative processes such as viscosity, resistivity and
anisotropic thermal conduction without operator splitting. Finally, we
illustrate an efficient treatment of pointlocal, potentially stiff source
terms over hierarchical nested grids by taking advantage of the adaptivity in
time. Several multidimensional benchmarks and applications to problems of
astrophysical relevance assess the potentiality of the AMR version of PLUTO in
resolving flow features separated by large spatial and temporal disparities.

We present and compare third as well as fifthorder accurate finite
difference schemes for the numerical solution of the compressible ideal MHD
equations in multiple spatial dimensions. The selected methods lean on four
different reconstruction techniques based on recently improved versions of the
weighted essentially nonoscillatory (WENO) schemes, monotonicity preserving
(MP) schemes as well as slopelimited polynomial reconstruction. The proposed
numerical methods are highly accurate in smooth regions of the flow, avoid loss
of accuracy in proximity of smooth extrema and provide sharp nonoscillatory
transitions at discontinuities. We suggest a numerical formulation based on a
cellcentered approach where all of the primary flow variables are discretized
at the zone center. The divergencefree condition is enforced by augmenting the
MHD equations with a generalized Lagrange multiplier yielding a mixed
hyperbolic/parabolic correction, as in Dedner et al. (J. Comput. Phys. 175
(2002) 645673). The resulting family of schemes is robust, costeffective and
straightforward to implement. Compared to previous existing approaches, it
completely avoids the CPU intensive workload associated with an elliptic
divergence cleaning step and the additional complexities required by staggered
mesh algorithms. Extensive numerical testing demonstrate the robustness and
reliability of the proposed framework for computations involving both smooth
and discontinuous features.

We investigate mode coupling in a two dimensional compressible disc with
radial stratification and differential rotation. We employ the global radial
scaling of linear perturbations and study the linear modes in the local
shearing sheet approximation. We employ a threemode formalism and study the
vorticity (W), entropy (S) and compressional (P) modes and their coupling
properties. The system exhibits asymmetric threemode coupling: these include
mutual coupling of S and Pmodes, S and Wmodes, and asymmetric coupling
between the W and Pmodes. Pmode perturbations are able to generate potential
vorticity through indirect threemode coupling. This process indicates that
compressional perturbations can lead to the development of vortical structures
and influence the dynamics of radially stratified hydrodynamic accretion and
protoplanetary discs.

Relativistic magnetized jets are key elements in Active Galactic Nuclei and
in other astrophysical environments. Their structure and evolution involves a
complex nonlinear physics that can be approached by numerical studies only.
Still, owing to a number of challenging computational aspects, only a few
numerical investigations have been undertaken so far. In this paper, we present
highresolution three dimensional numerical simulations of relativistic
magnetized jets carrying an initially toroidal magnetic field. The presence of
a substantial toroidal component of the field is nowadays commonly invoked and
held responsible for the process of jet acceleration and collimation. We find
that the typical nose cone structures, commonly observed in axisymmetric
twodimensional simulations, are not produced in the 3D case. Rather, the
toroidal field gives rise to strong current driven kink instabilities leading
to jet wiggling. However, it appears to be able to maintain an highly
relativistic spine along its full length. By comparing low and high resolution
simulations, we emphasize the impact of resolution on the jet dynamical
properties.

We present a fivewave Riemann solver for the equations of ideal relativistic
magnetohydrodynamics. Our solver can be regarded as a relativistic extension of
the fivewave HLLD Riemann solver initially developed by Miyoshi and Kusano for
the equations of ideal MHD. The solution to the Riemann problem is approximated
by a five wave pattern, comprised of two outermost fast shocks, two rotational
discontinuities and a contact surface in the middle. The proposed scheme is
considerably more elaborate than in the classical case since the normal
velocity is no longer constant across the rotational modes. Still, proper
closure to the RankineHugoniot jump conditions can be attained by solving a
nonlinear scalar equation in the total pressure variable which, for the chosen
configuration, has to be constant over the whole Riemann fan. The accuracy of
the new Riemann solver is validated against one dimensional tests and
multidimensional applications. It is shown that our new solver considerably
improves over the popular HLL solver or the recently proposed HLLC schemes.

We investigate the linear stability properties of the plane interface
separating two relativistic magnetized flows in relative motion. The two flows
are governed by the (special) relativistic equations for a magnetized perfect
gas in the infinite conductivity approximation. By adopting the vortexsheet
approximation, the relativistic magnetohydrodynamics equations are linearized
around the equilibrium state and the corresponding dispersion relation is
derived and discussed. The behavior of the configuration and the regimes of
instability are investigated following the effects of four physical parameters,
namely: the flow velocity, the relativistic and Alfv\'enic Mach numbers and the
inclination of the wave vector on the plane of the interface. From the
numerical solution of the dispersion relation, we find in general two separate
regions of instability, associated respectively with the slow and fast
magnetosonic modes. Modes parallel to the flow velocity are destabilized only
for sufficiently low magnetization. For the latter case, stabilization is
attained, additionally, at sufficiently large relativistic velocities between
the two flows in relative motion. The relevance of these results to the study
of the stability of astrophysical jets is briefly commented.

Strong observational evidence indicates that all extragalactic jets
associated with AGNs move at relativistic speed up to 100 pc  1 kpc scales
from the nucleus. At larger distances, reflecting the FanaroffRiley radio
source classification, we observe an abrupt deceleration in FRI jets while
relativistic motions persist up to Mpc scale in FRII. Moreover, VLBI
observations of some object like B2 1144+35, Mrk501 and M87 show limb
brightening of the jet radio emission at the parsec scale. This effect is
interpreted kinematically as due to the presence of a deboosted central spine
at high Lorentz factor and of a weakly relativistic external layer. In this
paper we investigate whether these effects can be interpreted by a breaking of
the collimated flow by external medium entrainment favored by shear
instabilities, namely KelvinHelmholtz instabilities. We examine in details the
physical conditions under which significant deceleration of a relativistic flow
is produced. We investigate the phenomenon by means of highresolution
threedimensional relativistic hydrodynamic simulations using the PLUTO code
for computational astrophysics. We find that the parameter of utmost importance
in determining the instability evolution and the entrainment properties is the
ambient/jet density contrast. We show that lighter jets suffer stronger slowing
down in the external layer than in the central part and conserve a central
spine at high Lorentz factor. Our model is verified by constructing synthetic
emission maps from the numerical simulations that compare reasonably well with
VLBI observations of the inner part of FRI sources.

Aims: We study the changes in the properties of turbulence driven by the
magnetorotational instability in a shearing box, as the computational domain
size in the radial direction is varied relative to the height Methods: We
perform 3D simulations in the shearing box approximation, with a net magnetic
flux, and we consider computational domains with different aspect ratios
Results: We find that in boxes of aspect ratio unity the transport of angular
momentum is strongly intermittent and dominated by channel solutions in
agreement with previous work. In contrast, in boxes with larger aspect ratio,
the channel solutions and the associated intermittent behavior disappear.
Conclusions: There is strong evidence that, as the aspect ratio becomes larger,
the characteristics of the solution become aspect ratio independent. We
conclude that shearing box calculations with aspect ratio unity or near unity
may introduce spurious effects.

We investigate the stability, nonlinear development and equilibrium structure
of vortices in a background shearing Keplerian flow. We make use of
highresolution global twodimensional compressible hydrodynamic simulations.
We introduce the concept of nonlinear adjustment to describe the transition of
unbalanced vortical fields to a longlived configuration. We discuss the
conditions under which vortical perturbations evolve into longlived persistent
structures and we describe the properties of these equilibrium vortices. The
properties of equilibrium vortices appear to be independent from the initial
conditions and depend only on the local disk parameters. In particular we find
that the ratio of the vortex size to the local disk scale height increases with
the decrease of the sound speed, reaching values well above the unity. The
process of spiral density wave generation by the vortex, discussed in our
previous work, appear to maintain its efficiency also at nonlinear amplitudes
and we observe the formation of spiral shocks attached to the vortex. The
shocks may have important consequences on the long term vortex evolution and
possibly on the global disk dynamics. Our study strengthens the arguments in
favor of anticyclonic vortices as the candidates for the promotion of planetary
formation. Hydrodynamic shocks that are an intrinsic property of persistent
vortices in compressible Keplerian flows are an important contributor to the
overall balance. These shocks support vortices against viscous dissipation by
generating local potential vorticity and should be responsible for the eventual
fate of the persistent anticyclonic vortices. Numerical codes have be able to
resolve shock waves to describe the vortex dynamics correctly.

We present a new numerical code, PLUTO, for the solution of hypersonic flows
in 1, 2 and 3 spatial dimensions and different systems of coordinates. The code
provides a multiphysics, multialgorithm modular environment particularly
oriented towards the treatment of astrophysical flows in presence of
discontinuities. Different hydrodynamic modules and algorithms may be
independently selected to properly describe Newtonian, relativistic, MHD or
relativistic MHD fluids. The modular structure exploits a general framework for
integrating a system of conservation laws, built on modern Godunovtype
shockcapturing schemes. Although a plethora of numerical methods has been
successfully developed over the past two decades, the vast majority shares a
common discretization recipe, involving three general steps: a piecewise
polynomial reconstruction followed by the solution of Riemann problems at zone
interfaces and a final evolution stage. We have checked and validated the code
against several benchmarks available in literature. Test problems in 1, 2 and 3
dimensions are discussed.