  Many astrophysical systems can only be accurately modelled when the behaviour
of their baryonic gas components is well understood. The residual distribution
(RD) family of partial differential equation (PDE) solvers produce approximate
solutions to the corresponding fluid equations. We present a new implementation
of the RD method. The solver efficiently calculates the evolution of the fluid,
with up to second order accuracy in both time and space, across an unstructured
triangulation, in both 2D and 3D. We implement a novel variable time stepping
routine, which applies a drifting mechanism to greatly improve the
computational efficiency of the method. We conduct extensive testing of the new
implementation, demonstrating its innate ability to resolve complex fluid
structures, even at very low resolution. We can resolve complex structures with
as few as 3-5 resolution elements, demonstrated by Kelvin-Helmholtz and Sedov
blast tests. We also note that we find cold cloud destruction time scales
consistent with those predicted by a typical PPE solver, albeit the exact
evolution shows small differences. The code includes three residual calculation
modes, the LDA, N and blended schemes, tailored for scenarios from smooth flows
(LDA), to extreme shocks (N), and both (blended). We compare our RD solver
results to state-of-the-art solvers used in other astrophysical codes,
demonstrating the competitiveness of the new approach, particularly at low
resolution. This is of particular interest in large scale astrophysical
simulations, where important structures, such as star forming gas clouds, are
often resolved by small numbers of fluid elements.
