  We present the first conservative finite volume numerical scheme for the
causal, stable relativistic Navier-Stokes equations developed by Bemfica,
Disconzi, Noronha, and Kovtun (BDNK). BDNK theory has arisen very recently as a
promising means of incorporating entropy-generating effects (viscosity, heat
conduction) into relativistic fluid models, appearing as a possible alternative
to the so-called M\"uller-Israel-Stewart (MIS) theory successfully used to
model quark-gluon plasma. The major difference between the two lies in the
structure of the system of PDEs: BDNK theory only has a set of conservation
laws, whereas MIS also includes a set of evolution equations for its
dissipative degrees of freedom. The simpler structure of the BDNK PDEs in this
respect allows for rigorous proofs of stability, causality, and hyperbolicity
in full generality which have as yet been impossible for MIS. To capitalize on
these advantages, we present the first fully conservative multi-dimensional
fluid solver for the BDNK equations suitable for physical applications. The
scheme includes a flux-conservative discretization, non-oscillatory
reconstruction, and a central-upwind numerical flux, and is designed to
smoothly transition to a high-resolution shock-capturing perfect fluid solver
in the inviscid limit. We assess the robustness of our new method in a series
of flat-spacetime tests for a conformal fluid, and provide a detailed
comparison with previous approaches of Pandya & Pretorius (2021).
