
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 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.

We investigate numerically Alfv\'en waves propagating along an axisymmetric
and nonisothermal solar flux tube embedded in the solar atmosphere. The tube
magnetic field is currentfree and diverges with height, and the waves are
excited by a periodic driver along the tube magnetic field lines. The main
results are that the two wave variables, the velocity and magnetic field
perturbations in the azimuthal direction, behave differently as a result of
gradients of physical parameters along the tube. To explain these differences
in the wave behavior, the time evolution of the wave variables and the
resulting cutoff period for each wave variable are calculated, and used to
determine regions in the solar chromosphere where strong wave reflection may
occur.

In the last decade, the relativistic magnetohydrodynamic (MHD) modelling of
pulsar wind nebulae, and of the Crab nebula in particular, has been highly
successful, with many of the observed dynamical and emission properties
reproduced down to the finest detail. Here, we critically discuss the results
of some of the most recent studies: namely the investigation of the origin of
the radio emitting particles and the quest for the acceleration sites of
particles of different energies along the termination shock, by using wisps
motion as a diagnostic tool; the study of the magnetic dissipation process in
high magnetization nebulae by means of new longterm threedimensional
simulations of the pulsar wind nebula evolution; the investigation of the
relativistic tearing instability in thinning current sheets, leading to fast
reconnection events that might be at the origin of the Crab nebula gammaray
flares.

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.

A significant fraction of OBtype, mainsequence massive stars are classified
as runaway and move supersonically through the interstellar medium (ISM). Their
strong stellar winds interact with their surroundings where the typical
strength of the local ISM magnetic field is about 3.57 microG, which can
result in the formation of bow shock nebulae. We investigate the effects of
such magnetic fields, aligned with the motion of the flow, on the formation and
emission properties of these circumstellar structures. Our axisymmetric,
magnetohydrodynamical simulations with opticallythin radiative cooling,
heating and anisotropic thermal conduction show that the presence of the
background ISM magnetic field affects the projected optical emission our bow
shocks at Ha and [OIII] lambda 5007 which become fainter by about 12 orders of
magnitude, respectively. Radiative transfer calculations against dust opacity
indicate that the magnetic field slightly diminishes their projected infrared
emission and that our bow shocks emit brightly at 60 micron. This may explain
why the bow shocks generated by ionizing runaway massive stars are often
difficult to identify. Finally, we discuss our results in the context of the
bow shock of Zeta Ophiuchi and we support the interpretation of its imperfect
morphology as an evidence of the presence of an ISM magnetic field not aligned
with the motion of its driving star.

We perform MHD modeling of a single bright coronal loop to include the
interaction with a nonuniform magnetic field. The field is stressed by random
footpoint rotation in the central region and its energy is dissipated into
heating by growing currents through anomalous magnetic diffusivity that
switches on in the corona above a current density threshold. We model an entire
single magnetic flux tube, in the solar atmosphere extending from the highbeta
chromosphere to the lowbeta corona through the steep transition region. The
magnetic field expands from the chromosphere to the corona. The maximum
resolution is ~30 km. We obtain an overall evolution typical of loop models and
realistic loop emission in the EUV and Xray bands. The plasma confined in the
flux tube is heated to active region temperatures (~3 MK) after ~2/3 hr.
Upflows from the chromosphere up to ~100 km/s fill the core of the flux tube to
densities above 10^9 cm^3. More heating is released in the low corona than the
high corona and is finely structured both in space and time.

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 present an interface between the (magneto) hydrodynamics code PLUTO and
the plasma simulation and spectral synthesis code CLOUDY. By combining these
codes, we constructed a new photoionization hydrodynamics solver: The
PLUTOCLOUDY Interface (TPCI), which is well suited to simulate
photoevaporative flows under strong irradiation. The code includes the
electromagnetic spectrum from Xrays to the radio range and solves the
photoionization and chemical network of the 30 lightest elements. TPCI follows
an iterative numerical scheme: First, the equilibrium state of the medium is
solved for a given radiation field by CLOUDY, resulting in a net radiative
heating or cooling. In the second step, the latter influences the (magneto)
hydrodynamic evolution calculated by PLUTO. Here, we validated the
onedimensional version of the code on the basis of four test problems:
Photoevaporation of a cool hydrogen cloud, cooling of coronal plasma, formation
of a Stroemgren sphere, and the evaporating atmosphere of a hot Jupiter. This
combination of an equilibrium photoionization solver with a general MHD code
provides an advanced simulation tool applicable to a variety of astrophysical
problems.

At least 5 per cent of the massive stars are moving supersonically through
the interstellar medium (ISM) and are expected to produce a stellar wind bow
shock. We explore how the mass loss and space velocity of massive runaway stars
affect the morphology of their bow shocks. We run twodimensional axisymmetric
hydrodynamical simulations following the evolution of the circumstellar medium
of these stars in the Galactic plane from the main sequence to the red
supergiant phase. We find that thermal conduction is an important process
governing the shape, size and structure of the bow shocks around hot stars, and
that they have an optical luminosity mainly produced by forbidden lines, e.g.
[OIII]. The Ha emission of the bow shocks around hot stars originates from near
their contact discontinuity. The H$\alpha$ emission of bow shocks around cool
stars originates from their forward shock, and is too faint to be observed for
the bow shocks that we simulate. The emission of opticallythin radiation
mainly comes from the shocked ISM material. All bow shock models are brighter
in the infrared, i.e. the infrared is the most appropriate waveband to search
for bow shocks. Our study suggests that the infrared emission comes from near
the contact discontinuity for bow shocks of hot stars and from the inner region
of shocked wind for bow shocks around cool stars. We predict that, in the
Galactic plane, the brightest, i.e. the most easily detectable bow shocks are
produced by highmass stars moving with small space velocities.

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.

Highorder reconstruction schemes for the solution of hyperbolic conservation
laws in orthogonal curvilinear coordinates are revised in the finite volume
approach. The formulation employs a piecewise polynomial approximation to the
zoneaverage values to reconstruct left and right interface states from within
a computational zone to arbitrary order of accuracy by inverting a
Vandermondelike linear system of equations with spatially varying
coefficients. The approach is general and can be used on uniform and
nonuniform meshes although explicit expressions are derived for polynomials
from second to fifth degree in cylindrical and spherical geometries with
uniform grid spacing. It is shown that, in regions of large curvature, the
resulting expressions differ considerably from their Cartesian counterparts and
that the lack of such corrections can severely degrade the accuracy of the
solution close to the coordinate origin. Limiting techniques and monotonicity
constraints are revised for conventional reconstruction schemes, namely, the
piecewise linear method (PLM), thirdorder weighted essentially nonoscillatory
(WENO) scheme and the piecewise parabolic method (PPM).
The performance of the improved reconstruction schemes is investigated in a
number of selected numerical benchmarks involving the solution of both scalar
and systems of nonlinear equations (such as the equations of gas dynamics and
magnetohydrodynamics) in cylindrical and spherical geometries in one and two
dimensions. Results confirm that the proposed approach yields considerably
smaller errors, higher convergence rates and it avoid spurious numerical
effects at a symmetry axis.

The expansion of coronal loops in the transition region may considerably
influence the diagnostics of the plasma emission measure. The cross sectional
area of the loops is expected to depend on the temperature and pressure, and
might be sensitive to the heating rate. The approach here is to study the area
response to slow changes in the coronal heating rate, and check the current
interpretation in terms of steady heating models. We study the area response
with a timedependent 2D MHD loop model, including the description of the
expanding magnetic field, coronal heating and losses by thermal conduction and
radiation from optically thin plasma. We run a simulation for a loop 50 Mm long
and quasistatically heated to about 4 MK. We find that the area can change
substantially with the quasisteady heating rate, e.g. by ~40% at 0.5 MK as the
loop temperature varies between 1 and 4 MK, and, therefore, affects the
interpretation of DEM(T) curves.

According to the magnetospheric accretion scenario, young lowmass stars are
surrounded by circumstellar disks which they interact with through accretion of
mass. The accretion builds up the star to its final mass and is also believed
to power the mass outflows, which may in turn have a significant role in
removing the excess angular momentum from the stardisk system. Although the
process of mass accretion is a critical aspect of star formation, some of its
mechanisms are still to be fully understood. On the other hand, strong flaring
activity is a common feature of young stellar objects (YSOs). In the Sun, such
events give rise to perturbations of the interplanetary medium. Similar but
more energetic phenomena occur in YSOs and may influence the circumstellar
environment. In fact, a recent study has shown that an intense flaring activity
close to the disk may strongly perturb the stability of circumstellar disks,
thus inducing mass accretion episodes (Orlando et al. 2011). Here we review the
main results obtained in the field and the future perspectives.

Astronomical observations, analytical solutions and numerical simulations
have provided the building blocks to formulate the current theory of young
stellar object jets. Although each approach has made great progress
independently, it is only during the last decade that significant efforts are
being made to bring the separate pieces together. Building on previous work
that combined analytical solutions and numerical simulations, we apply a
sophisticated cooling function to incorporate optically thin energy losses in
the dynamics. On the one hand, this allows a selfconsistent treatment of the
jet evolution and on the other, it provides the necessary data to generate
synthetic emission maps. Firstly, analytical disk and stellar outflow solutions
are properly combined to initialize numerical twocomponent jet models inside
the computational box. Secondly, magnetohydrodynamical simulations are
performed in 2.5D, following properly the ionization and recombination of a
maximum of $29$ ions. Finally, the outputs are postprocessed to produce
artificial observational data. The first twocomponent jet simulations, based
on analytical models, that include ionization and optically thin radiation
losses demonstrate promising results for modeling specific young stellar object
outflows. The generation of synthetic emission maps provides the link to
observations, as well as the necessary feedback for the further improvement of
the available models.

We investigate the dynamical propagation of the SouthEast jet from the Crab
pulsar interacting with supernova ejecta by means of threedimensional
relativistic MHD numerical simulations with the PLUTO code.
The initial jet structure is set up from the inner regions of the Crab
Nebula.
We study the evolution of hot, relativistic hollow outflows initially
carrying a purely azimuthal magnetic field.
Our jet models are characterized by different choices of the outflow
magnetization ($\sigma$ parameter) and the bulk Lorentz factor ($\gamma_{j}$).
We show that the jet is heavily affected by the growth of currentdriven kink
instabilities causing considerable deflection throughout its propagation
length.
This behavior is partially stabilized by the combined action of larger flow
velocities and/or reduced magnetic field strengths.
We find that our best jet models are characterized by relatively large values
of $\sigma$ ($\gtrsim 1$) and small values of $\gamma_{j}\simeq 2$.
Our results are in good agreement with the recent Xray (\textit{Chandra})
data of the Crab Nebula SouthEast jet indicating that the jet changes
direction of propagation on a time scale of the order of few years.
The 3D models presented here may have important implications in the
investigation of particle acceleration in relativistic outflows.

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.

It is a well established fact that some YSO jets (e.g. RW Aur) display
different propagation speeds between their blue and red shifted parts, a
feature possibly associated with the central engine or the environment in which
the jet propagates. In order to understand the origin of asymmetric YSO jet
velocities, we investigate the efficiency of two candidate mechanisms, one
based on the intrinsic properties of the system and one based on the role of
the external medium. In particular, a parallel or antiparallel configuration
between the protostellar magnetosphere and the disk magnetic field is
considered and the resulting dynamics are examined both in an ideal and a
resistive magnetohydrodynamical (MHD) regime. Moreover, we explore the effects
of a potential difference in the pressure of the environment, as a consequence
of the nonuniform density distribution of molecular clouds. Ideal and
resistive axisymmetric numerical simulations are carried out for a variety of
models, all of which are based on a combination of two analytical solutions, a
disk wind and a stellar outflow. We find that jet velocity asymmetries can
indeed occur both when multipolar magnetic moments are present in the stardisk
system as well as when nonuniform environments are considered. The latter case
is an external mechanism that can easily explain the large time scale of the
phenomenon, whereas the former one naturally relates it to the YSO intrinsic
properties. [abridged]

Explicit numerical computations of superfast differentially rotating disks
are subject to the timestep constraint imposed by the Courant condition. When
the bulk orbital velocity largely exceeds any other wave speed the time step is
considerably reduced and a large number of steps may be necessary to complete
the computation.
We present a robust numerical scheme to overcome the Courant limitation by
extending the algorithm previously known as FARGO (Fast Advection in Rotating
Gaseous Objects) to the equations of magnetohydrodynamics (MHD). The proposed
scheme conserves total angular momentum and energy to machine precision and
works in Cartesian, cylindrical, or spherical coordinates. The algorithm is
implemented in the PLUTO code for astrophysical gasdynamics and is suitable for
local or global simulations of accretion or protoplanetary disk models.
By decomposing the total velocity into an average azimuthal contribution and
a residual term, the algorithm solves the MHD equations through a linear
transport step in the orbital direction and a standard nonlinear solver applied
to the MHD equations written in terms of the residual velocity. Since the
former step is not subject to any stability restriction, the Courant condition
is computed only in terms of the residual velocity, leading to substantially
larger time steps. The magnetic field is advanced in time using the constrained
transport method in order to preserve the divergencefree condition.
Conservation of total energy and angular momentum is enforced at the discrete
level by properly expressing the source terms in terms of upwind fluxes
available during the standard solver.
Our results show that applications of the proposed orbitaladvection scheme
to problems of astrophysical relevance provides, at reduced numerical cost,
equally accurate and less dissipative results than standard timemarching
schemes.

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.