  The relativistic hydrodynamics (RHD) equations have three crucial intrinsic
physical constraints on the primitive variables: positivity of pressure and
density, and subluminal fluid velocity. However, numerical simulations can
violate these constraints, leading to nonphysical results or even simulation
failure. Designing genuinely physical-constraint-preserving (PCP) schemes is
very difficult, as the primitive variables cannot be explicitly reformulated
using conservative variables due to relativistic effects. In this paper, we
propose three efficient Newton--Raphson (NR) methods for robustly recovering
primitive variables from conservative variables. Importantly, we rigorously
prove that these NR methods are always convergent and PCP, meaning they
preserve the physical constraints throughout the NR iterations. The discovery
of these robust NR methods and their PCP convergence analyses are highly
nontrivial and technical. As an application, we apply the proposed NR methods
to design PCP finite volume Hermite weighted essentially non-oscillatory
(HWENO) schemes for solving the RHD equations. Our PCP HWENO schemes
incorporate high-order HWENO reconstruction, a PCP limiter, and
strong-stability-preserving time discretization. We rigorously prove the PCP
property of the fully discrete schemes using convex decomposition techniques.
Moreover, we suggest the characteristic decomposition with rescaled
eigenvectors and scale-invariant nonlinear weights to enhance the performance
of the HWENO schemes in simulating large-scale RHD problems. Several demanding
numerical tests are conducted to demonstrate the robustness, accuracy, and high
resolution of the proposed PCP HWENO schemes and to validate the efficiency of
our NR methods.
