For simplicity I'll assume $u$ is scalar valued, but I am pretty sure the discussion also works for $u$ that is a section of some vector bundle over $M$ (if the wave operator is quasidiagonal). Additionally, the computations are all essentially local, so I'll forget that $M$ is a manifold and just work in local coordinates.
We have a domain $\Omega\subseteq \mathbb{R}^n$, and are given a matrix-valued function $\mathbf{g}: \Omega\times \mathbb{R}\to \mathrm{Lor}(n)$, where $\mathrm{Lor}(n)$ denote the $n\times n$ symmetric matrices with Lorentzian signature. We are also given a function $F: \Omega \times \mathbb{R}\times \mathbb{R}^n \to \mathbb{R}$, and we wish to study the partial differential equation
$$ \frac{1}{\sqrt{|\mathbf{g}(x,u)|}} \partial_i \Big[(\mathbf{g}(x,u)^{-1})^{ij} \sqrt{|\mathbf{g}(x,u)|} \partial_j u \Big] = F(x,u, \nabla u) \tag{*}\label{eq:basic}$$
The division by the scalar $\sqrt{|\mathbf{g}|}$ is inconsequential; we can simply set $\tilde{F} = F\cdot \sqrt{|\mathbf{g}|}$. Our first observation is that if $u\in C^2(\Omega)$ is a classical solution to \eqref{eq:basic}, then for any $\phi \in C^1_c(\Omega)$ we have that
$$ \int_{\Omega} (\mathbf{g}^{-1})^{ij} \sqrt{|\mathbf{g}|} \partial_i \phi \partial_j u ~dx + \int_{\Omega} \phi \tilde{F} ~dx = 0 \tag{**}\label{eq:weak} $$
We can regard \eqref{eq:weak} as the "weak formulation" of \eqref{eq:basic}: provided we have $u$ defined in such a way that $\tilde{F}$ makes sense as an element in $(C^1_c(\Omega))'$ and that $(\mathbf{g}^{-1})^{ij} \sqrt{|\mathbf{g}|} \partial_i u$ makes sense as an element of $(C^0_c(\Omega))'$, we can say that "$u$ is a weak solution iff \eqref{eq:weak} is satisfied".
One can discuss this in $BV$-based spaces, but one can also derive the Rankine-Hugoniot conditions without thinking about solutions with such low regularity.
(That $BV$-based spaces is often used when studying conservation laws is because when we are in $n = 2$ (so spatial dimension 1) there is a robust theory of existence uniqueness in this class. But since such a theory doesn't exist in higher dimensions, I won't bother here.)
Since we are interested in shocks, let's make (the fairly strong) assumption that the shock surface is a regular hypersurface $\Sigma$ of $\Omega$ that divides $\Omega$ into two open subsets $\Omega_\pm$. Suppose we have $u:C^2(\Omega_+\cup \Omega_-)$, such that $u|_{\Omega_+}$ has a continuous extension to $\Omega_+\cap \Sigma$, and similarly $\partial u|_{\Omega_+}$; ditto the restrictions to $\Omega_-$. We will denote by $u_+:\Omega_+\cap \Sigma\to \mathbb{R}$ said extension and similarly $u_-:\Omega_-\cap \Sigma\to \mathbb{R}$.
Now let $n$ be a normal co-vector field for $\Sigma$. Then it is not too hard to see that
If $u$ is a function as above, then $u$ solves \eqref{eq:weak} if and only if $u|_{\Omega_\pm}$ solves \eqref{eq:basic} and that at every $p\in \Sigma$, we have $$(\mathbf{g}^{-1}(p,u_+(p)))^{ij} \sqrt{|\mathbf{g}(p,u_+(p))|} n_i(p) \nabla_j u_+(p) = (\mathbf{g}^{-1}(p,u_-(p)))^{ij} \sqrt{|\mathbf{g}(p,u_-(p))|} n_i(p) \nabla_j u_-(p) \tag{RH} \label{eq:RH}$$
Sketch of proof: if $\mathrm{supp}(\phi)$ is disjoint from $\Sigma$, integrate by parts. This prove the requirement of the first condition. If $\mathrm{supp}(\phi)$ intersects $\Sigma$, using that $\Sigma$ has measure zero replace \eqref{eq:weak} with integration over $\Omega\setminus \Sigma$. Our conditions on $u$ allows us to integrate by parts separately on $\Omega_+$ and $\Omega_-$; since \eqref{eq:basic} holds away from $\Sigma$ we just need the boundary terms to be equal.
Note: The equation \eqref{eq:RH} is the Rankine-Hugoniot condition. This proof is fundamentally identical to the proof of the Rankine-Hugoniot condition for $BV_{\mathrm{loc}}$ solutions to the first order hyperbolic systems; the $BV$ case just has more technical definitions.
The entropy condition is something different.
One of the uses of the entropy condition is to guarantee well-posedness of the initial value problem; it rules out the spontaneous formation of rarefaction waves. This can be regarded as a causality requirement.
In terms of the causality requirement, this can be generalized to the wave equations case pretty easily. For this, however, we need an "arrow of time".
Assumption (this happens to hold for irrotational homentropic fluid systems too) there exists a (non-vanishing) vector field $T$ on $M$ such that for every metric in the range of $\mathbf{g}$, the vector field $T$ is a time-like vector field.
Then a reasonable causality condition is this: the admissible shock surfaces $\Sigma$ are such that $T$ is transverse to $\Sigma$. Furthermore, assuming $T$ points from $\Omega_-$ to $\Omega_+$ across $\Sigma$, then $\Sigma$ is space-like relative to $\mathbf{g}(x,u_-(x))$, and $\Sigma$ is time-like relative to $\mathbf{g}(x,u_+(x))$.
(Basically, the idea is that the shock front should exceed the speed of sound of the fluid "in front of it", hence it should be a space-like surface with respect to that fluid. On the other hand, if the shock front were also space-like relative to the fluid that is behind it, then the initial value problem is well-posed coming from it. This would necessarily imply the non-uniqueness of the shock front. So we should assume the opposite condition.)
You can see that this also agrees with the Lax conditions for scalar conservation laws in 1D.
Off hand I don't remember if additional entropy-type conditions are needed beyond this.
