Conservative finite volume scheme for first-order viscous relativistic hydrodynamics

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üller-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 partial differential equations (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 multidimensional fluid solver for the BDNK equations suitable for physical applications. The scheme includes a flux-conservative discretization, nonoscillatory 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 and Pretorius [Phys. Rev. D 104, 023015 (2021)].

Figure 1

Discrete conservation of $T^{tt}$ across the spatial domain for a simulation starting from Gaussian initial data (48) with ${\eta }_0=0.2$ for the finite volume scheme presented here (FV) as well as the semi-finite-difference scheme of [36] (FD). As expected, the finite volume scheme conserves $T^{tt}$ up to machine precision $\sim {10}^{-15}$, until the fluid pulse reaches the boundary at the time marked by the light red dotted line. The semi-finite-difference scheme of [36] conserves $T^{tt}$ only up to the level of truncation error, which in this case is $\sim {10}^{-3}$.

Figure 2

Density ($n$) evolution of viscous rotor initial data (49) as a function of time (columns) for three different viscosities: ${\eta }_0=0$, 0.01, 0.2 in rows, from top to bottom. In the inviscid simulation, the cylinder of fluid is Kelvin-Helmholtz unstable and forms vortices which are not present in the viscous cases. At intermediate viscosity (middle row), the fluid experiences a shearing force which distorts the bar of overdensity present in the initial data, before the cylinder stops rotating entirely around $t\sim 12$. At the highest viscosity shown (bottom), the cylinder rotates only about 20° before stopping at $t\sim 5$.

Figure 3

Integral of the absolute value of the constraint (42) over the domain for viscous rotor initial data with ${\eta }_0=0.2$. The value of ${\epsilon }_W$ determines how strongly the smoothness of the candidate ENO stencils impacts the nonlinear weights; small values of ${\epsilon }_W$ imply strong sensitivity to nonsmoothness, and large values imply insensitivity (and as a result give a fixed fourth-order derivative stencil). Hence, for smaller ${\epsilon }_W$ one finds larger violations of the constraint (42), which converge away with resolution (the solid lines range over $N_x=2^7,2^8,2^9$, with lighter shades representing higher resolutions). In the ${\epsilon }_W\to \infty$ limit, constraint violation approaches machine precision (cf. the ${\epsilon }_W={10}^{15}$ case).

Figure 4

Comparison of solutions for $\epsilon$ starting from shock tube initial data (51) at three successive resolutions for the semi-finite-difference scheme of [36] (left, FD) versus the finite volume scheme presented here (right, FV) at $t\sim 43$ for ${\eta }_0=0.2$. The FD scheme has oscillations near the shock front which quickly converge away with resolution, as well as grid-scale “sawtooth” oscillations that developed early on near the origin (the $t=0$ location of the shock front) and do not converge away as rapidly with resolution. These features do not appear in the figures of [36] because the discontinuities there are smaller in amplitude, leading to oscillations small enough to be tamed by applying Kreiss-Oliger dissipation; said dissipation is not strong enough to remove the oscillations for the initial data (51), and we choose not to apply artificial dissipation in either scheme throughout this work. The FV solutions are free of noticeable oscillations, and the $N_x=2^9,2^{10},2^{11}$ curves all overlap at the resolution of the plot.

Figure 5

Solution for the log of the energy density $\xi$ for the 2D oblique shock wave initial data (53) at $t\sim 220$ for ${\eta }_0=0.2$. Note that the solution is nonoscillatory, even though there is an order unity jump in $\xi$ (corresponding to a jump of $\sim 50$ in $\epsilon$) which is not aligned with the numerical grid.

Figure 6

Illustration of the algorithm used to capture the perfect fluid limit for steady-state shock wave initial data (56) at ${\eta }_0=0.2$. The dashed line (top) is the solution for this set of initial data at late times, constructed using the BDNK primitive variable solution (39) everywhere. When the adaptive primitive variable solver is used, (39) is only used in the gray region, where the shade of gray corresponds to the value of the viscous tolerance ${\mathrm{\Delta }}_{\eta }$ shown in the legend. For large values of this tolerance, only regions with very steep gradients are identified as being nonequilibrium, and (39) is only used in a small sliver of the solution (and the perfect fluid primitive variable solution (29) is used elsewhere). This induces significant errors (bottom) when compared to the solution where only (39) is used. Shrinking the viscous tolerance ${\mathrm{\Delta }}_{\eta }$ results in more of the nonequilibrium region being identified as such by the algorithm and gives successively smaller errors when compared to the dashed (BDNK-only) solution. For ${\mathrm{\Delta }}_{\eta }\lesssim {10}^{-7}$, the error drops to machine precision.

Figure 7

Evolution of Kelvin-Helmholtz-unstable initial data (57) for the density $n$ at three different viscosities in columns, from left to right: ${\eta }_0=0,4\times {10}^{-4},{10}^{-3}$, at $t=11$ (top row) and $t=31$ (bottom row). Viscosity has a clear effect on both the early- and late-time state of the fluid; at $t=11$ it determines the amount of growth of the perturbation of low-density fluid (dark blue) into the high-density (yellow) region. For the two lower-viscosity cases (left two columns), long-lived vortices form out of these perturbations. At high viscosity, no clear vortex has formed, instead the perturbation has been sheared into a long, thin mixed layer.

Figure 8

Snapshot of the density $n$ at $t=11$ from the Kelvin-Helmholtz simulation for ${\eta }_0=4\times {10}^{-4}$, as a function of resolution. For the two lower-resolution panels (top row), the physical viscosity is small enough that the BDNK primitive variable solution is numerically unstable, and we use the adaptive primitive variable algorithm with tolerances ${\mathrm{\Delta }}_{\eta }={10}^{-3},{10}^{-4}$, respectively. At these tolerances, the perfect fluid primitive variable solution is used across the entire grid for most of the simulation after the first few time steps; despite this, the solution still converges to the correct viscous solution and is noticeably different from the inviscid solution (top left of Fig. 7) because most of the dissipation comes from the viscous fluxes rather than the primitive variable solution.

Figure 9

Convergence plots corresponding to an independent (Crank-Nicolson second-order finite-difference) discretization of the $x$ component of (1), ${\nabla }_c{T}^{cx}=0$, for the 2D simulations shown above. Leftmost column: plots of $Q_N(t)$ for the viscous rotor simulations showing that the curves approach second-order convergence as resolution increases (shown with finer grids in successively darker colors, corresponding to resolutions $N_x=2^8,2^9,2^{10}$). Middle column: similar plots of $Q_N(t)$ for the Kelvin-Helmholtz-unstable initial data. Rightmost panel: plots of the residual ${\nabla }_c{T}^{cx}$ for the 2D oblique shock initial data, at the time shown in Fig. 5, with the same color coding by resolution, scaled such that all curves should overlap if they are converging at second order, e.g., the $N_x=2^9$ curve is multiplied by 4 and the $N_x=2^{10}$ curve is multiplied by 16. The three panels show successive slices through the domain at constant $x$, and the top two show convergence at the expected order (all curves overlap). The bottom curve shows a slice through the shock wave and converges roughly at the expected order everywhere outside the spikes which appear at the shock fronts. Increasing resolution should produce taller, thinner spikes at the shock wave until it is finally resolved and the solution begins converging there at second order.