  We present an extension of the Piecewise Parabolic Method to special
relativistic fluid dynamics in multidimensions. The scheme is conservative,
dimensionally unsplit, and suitable for a general equation of state. Temporal
evolution is second-order accurate and employs characteristic projection
operators; spatial interpolation is piece-wise parabolic making the scheme
third-order accurate in smooth regions of the flow away from discontinuities.
The algorithm is written for a general system of orthogonal curvilinear
coordinates and can be used for computations in non-cartesian geometries. A
non-linear iterative Riemann solver based on the two-shock approximation is
used in flux calculation. In this approximation, an initial discontinuity
decays into a set of discontinuous waves only implying that, in particular,
rarefaction waves are treated as flow discontinuities. We also present a new
and simple equation of state which approximates the exact result for the
relativistic perfect gas with high accuracy. The strength of the new method is
demonstrated in a series of numerical tests and more complex simulations in
one, two and three dimensions.
