  Cosmic rays (CRs) are frequently modeled as an additional fluid in
hydrodynamic (HD) and magnetohydrodynamic (MHD) simulations of astrophysical
flows. The standard CR two-fluid model is described in terms of three
conservation laws (expressing conservation of mass, momentum and total energy)
and one additional equation (for the CR pressure) that cannot be cast in a
satisfactory conservative form. The presence of non-conservative terms with
spatial derivatives in the model equations prevents a unique weak solution
behind a shock. We investigate a number of methods for the numerical solution
of the two-fluid equations and find that, in the presence of shock waves, the
results generally depend on the numerical details (spatial reconstruction, time
stepping, the CFL number, and the adopted discretization). All methods converge
to a unique result if the energy partition between the thermal and non-thermal
fluids at the shock is prescribed using a subgrid prescription. This highlights
the non-uniqueness problem of the two-fluid equations at shocks. From our
numerical investigations, we report a robust method for which the solutions are
insensitive to the numerical details even in absence of a subgrid prescription,
although we recommend a subgrid closure at shocks using results from kinetic
theory. The subgrid closure is crucial for a reliable post-shock solution and
also its impact on large scale flows because the shock microphysics that
determines CR acceleration is not accurately captured in a fluid approximation.
Critical test problems, limitations of fluid modeling, and future directions
are discussed.
