  We discuss new methods to integrate the cosmic ray (CR) evolution equations
coupled to magneto-hydrodynamics (MHD) on an unstructured moving mesh, as
realised in the massively parallel AREPO code for cosmological simulations. We
account for diffusive shock acceleration of CRs at resolved shocks and at
supernova remnants in the interstellar medium (ISM), and follow the advective
CR transport within the magnetised plasma, as well as anisotropic diffusive
transport of CRs along the local magnetic field. CR losses are included in
terms of Coulomb and hadronic interactions with the thermal plasma. We
demonstrate the accuracy of our formalism for CR acceleration at shocks through
simulations of plane-parallel shock tubes that are compared to newly derived
exact solutions of the Riemann shock tube problem with CR acceleration. We find
that the increased compressibility of the post-shock plasma due to the produced
CRs decreases the shock speed. However, CR acceleration at spherically
expanding blast waves does not significantly break the self-similarity of the
Sedov-Taylor solution; the resulting modifications can be approximated by a
suitably adjusted, but constant adiabatic index. In first applications of the
new CR formalism to simulations of isolated galaxies and cosmic structure
formation, we find that CRs add an important pressure component to the ISM that
increases the vertical scale height of disk galaxies, and thus reduces the star
formation rate. Strong external structure formation shocks inject CRs into the
gas, but the relative pressure of this component decreases towards halo centres
as adiabatic compression favours the thermal over the CR pressure.
