SRHD

Derivation:
for physical derivation:
2008 Font "Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity"

relativisitc metric
$g_{\mu\nu}$ = metric tensor.
ADM metric:
lapse: $ \alpha $
shift spatial components: $ \beta^i $
spatial metric spatial components: $ \gamma_{ij} $
$||g_{\mu\nu}|| = \overset{\mu(i)\downarrow\nu(j)\rightarrow}{\left[ \matrix{-\alpha^2 + \beta_k \beta^k & \beta_j \\ \beta_i & \gamma_{ij}} \right]}$
$||g^{\mu\nu}|| = \overset{\mu(i)\downarrow\nu(j)\rightarrow}{\left[ \matrix{-1/\alpha^2 & \beta^j / \alpha^2 \\ \beta^i / \alpha^2 & \gamma^{ij} - \beta^i \beta^j / \alpha^2} \right]}$
where $\beta_j = \gamma_{ij} \beta^i = \gamma_{\mu j} \beta^\mu = g_{\mu j} \beta^j $ courtesy of the fact that $ g_{ij} = \gamma_{ij} $ and $ \beta^t = 0 $.
ADM metric determinant: $ g = \alpha \sqrt{\gamma} $
ADM spatial metric determinant: $ \gamma = det||\gamma_{ij}||$
normal in the ADM metric:
$n_\mu = \overset{\mu(j)\rightarrow}{\left[ \matrix{ -\alpha & 0 } \right]} $
$n^\mu = \overset{\mu(i)\downarrow}{\left[ \matrix{ 1/\alpha \\ -\beta^i/\alpha } \right]} $
For SR the metric is set to $g_{\mu\nu} = \eta_{\mu\nu} = diag(-1,1,1,1)$
so $\alpha = 1$ and $\beta^i = 0$ and $\gamma_{ij} = \delta_{ij} = diag(1,1,1)$

some more variables
$\rho$ = rest-mass density
$u^\mu$ = 4-velocity
$v^i$ = Newtonian velocity
$\epsilon$ = specific internal energy

current of rest-mass
$J^\mu = \rho u^\mu$ = current of rest-mass
$\nabla_\mu J^\mu = 0$ shows current of rest-mass is conserved

pressure
$p$ = rest-mass pressure
equation of state:
$p = p(\rho, \epsilon)$
pressure for ideal gas:
$\gamma$ = heat capacity ratio
$p = (\gamma - 1) \rho \epsilon$

enthalpy
In non-relativistic case:
$h = \epsilon + p / \rho$
this page (which doesn't cite its sources) says:
$\epsilon$ = the rest-mass specific internal energy
$\rho \epsilon_R = \rho c^2 + \rho \epsilon$, for $\epsilon_R$ = the relativistic specific internal energy
In the relativistic case:
$h = \epsilon_R + p / \rho = c^2 + \epsilon + p / \rho$
Note that this source says the speed of sound $c_s = \frac{\partial p}{\partial e_R}$ while my SRHD sources say $c_s = \frac{\partial p}{\partial e}$
rest-mass internal specific enthalpy:
$h = 1 + \epsilon + p / \rho$
$h = 1 + \gamma \epsilon$

stress-energy tensor
$T^{\mu\nu}$ = stress-energy / energy-momentum tensor
$\nabla_\mu T^{\mu\nu} = 0$ shows stress-energy is conserved
$T^{\mu\nu} = \rho h u^\mu u^\nu + p g^{\mu\nu}$ = definition of stress-energy tensor for perfect fluids
(This is also given in plenty of other sources, M.T.W., Alcubierre, Baumgarte & Shapiro, etc.)

speed of sound
$\chi = \partial p / \partial \rho |_\epsilon$
$\kappa = \partial p / \partial \epsilon |_\rho$
$h c_s^2 = \chi + p/\rho^2 \kappa$

relativistic Lorentz factor
Let $v^i$ be our spatial velocity, related to $u^i$ by $v^i = \frac{u^i}{\alpha u^t} + \frac{\beta^i}{\alpha}$
Let $W = \alpha u^t = 1/\sqrt{1 - v^i v^j \gamma_{ij}}$

For SRHD:
$W = 1 / \sqrt{1 - v^2}$
$W^2 = 1 / (1 - v^2)$
$1 / W^2 = 1 - v^2$
$v^2 = 1 - 1 / W^2$

primitive variables
$W_I = \downarrow I \left[\matrix{\rho \\ v_i \\ \epsilon}\right]$ = primitive variables

conservative variables
These variables appear just as they do in the non-relativistic Euler fluid equations, however $\rho$ is now densitized by $W$, and $S^i$ and $\tau$ are now densitized by an additional $h W$
$D = -J^\mu n_\mu = \alpha J^0 = \alpha \rho u^0 = \rho W$ = relativistic density
$S^i = -T^{\mu i} n_\mu = \alpha T^{0i} = \alpha (\rho h u^0 u^i + p g^{0i}) = \rho h W^2 v^i $ = relativistic momentum
$E = T^{\mu\nu} n_\mu n_\nu = \alpha^2 T^{00} = \alpha^2 (\rho h (u^0)^2 + p g^{00}) = \rho h W^2 - p$ = relativistic total energy density
$\tau = E - D = \rho h W^2 - p - D$ = relativistic total energy density minus rest-mass density

$U_I = \downarrow I \left[\matrix{D \\ S_i \\ \tau }\right]$ = conservative variables

$F_{Ij} = \downarrow I \left[\matrix{D v_j \\ S_i v_j - \delta_{ij} p \\ S_j - D v_j}\right]$ = flux vector

$\partial_\rho$:
$P_{,\rho} = (\gamma - 1) \epsilon$
$h_{,\rho} = P_{,\rho} / \rho - p / \rho^2 = ((\gamma - 1) \rho \epsilon - p) / \rho^2 = 0$
$D_{,\rho} = W$
$S_{i,\rho} = h W^2 v_i$
$\tau_{,\rho} = h W^2 - P_{,\rho} - D_{,\rho} = h W^2 - (\gamma - 1) \epsilon - W = W (h W - 1) - (\gamma - 1) \epsilon$

$\partial_{v_k}$:
$W_{,v_k} = \partial_{v_k} (\sqrt{1-v^2})^{-1} = -\frac{1}{2} (\sqrt{1-v^2})^{-3} \cdot -2 v^k = \frac{v^k}{(\sqrt{1-v^2})^3} = v^k W^3$
$D_{,v_k} = \rho W_{,v_k} = \rho W^3 v^k = W S^k / h$
$S_{i,v_k} = \rho h (2 W W_{,v_k} v_i + W^2 \delta_i^k) = \rho h W^2 (2 W^2 v_i v^k + \delta_i^k)$
$\tau_{,v_k} = \rho h 2 W W_{,v_k} - D_{,v_k} = 2 \rho h W^4 v^k - \rho W^3 v^k = \rho W^3 v^k (2 h W - 1)$

$\partial_{,\epsilon}$:
$P_{,\epsilon} = (\gamma - 1) \rho$
$h_{,\epsilon} = 1 + P_{,\epsilon} / \rho = 1 + \gamma - 1 = \gamma$
$S_{i,\epsilon} = \rho W^2 v_i h_{,\epsilon} = \gamma \rho W^2 v_i = \gamma S_i / h$
$\tau_{,\epsilon} = \rho h_{,\epsilon} W^2 - P_{,\epsilon} = \gamma \rho W^2 - (\gamma - 1) \rho = \rho (\gamma (W^2 - 1) + 1)$

change in conserved variables wrt primitive variables:
$\frac{\partial U_I}{\partial W_J} = \downarrow I \overset{ \rightarrow J }{ \left[\matrix{ D_{,\rho} & D_{,v_j} & D_{,\epsilon} \\ S_{i,\rho} & S_{i,v_j} & S_{i,\epsilon} \\ \tau_{,\rho} & \tau_{,v_j} & \tau_{,\epsilon} }\right] }$

$\frac{\partial U_I}{\partial W_J} = \downarrow I \overset{ \rightarrow J }{ \left[\matrix{ W & \rho W^3 v_j & 0 \\ h W^2 v_i & \rho h W^2 (2 W^2 v_i v_j + \delta_{ij}) & \gamma \rho W^2 v_i \\ W (h W - 1) - (\gamma - 1) \epsilon & \rho W^3 v_j (2 h W - 1) & \rho (\gamma W^2 - (\gamma - 1)) }\right] }$

change in primitive variables wrt conserved variables:

$\frac{\partial W_J}{\partial U_I} = [\frac{\partial U_I}{\partial W_J}]^{-1}$

Inverting with Gauss-Jordan elimination:
$\left[\begin{matrix} W & \rho W^3 v_j & 0 & | & 1 & 0 & 0 \\ h W^2 v_i & \rho h W^2 (2 W^2 v_i v_j + \delta_{ij}) & \gamma \rho W^2 v_i & | & 0 & \delta_{ij} & 0 \\ W (h W - 1) - (\gamma - 1) \epsilon & \rho W^3 v_j (2 h W - 1) & \rho (\gamma W^2 - (\gamma - 1)) & | & 0 & 0 & 1 \end{matrix}\right]$
$A \rightarrow \frac{1}{W} A$
$B \rightarrow B - h W v_i A$
$C \rightarrow C - A ( h W - 1 - \frac{1}{W} (\gamma - 1) \epsilon)$
$\left[\begin{matrix} 1 & \rho W^2 v_j & 0 & | & \frac{1}{W} & 0 & 0 \\ 0 & \rho h W^2 (2 W^2 v_i v_j + \delta_{ij}) - \rho h W^4 v_i v_j & \gamma \rho W^2 v_i & | & - h W v_i & \delta_{ij} & 0 \\ 0 & \rho W^3 v_j (2 h W - 1) - \rho W^3 v_j ( h W - 1 - \frac{1}{W} (\gamma - 1) \epsilon) & \rho (\gamma W^2 - (\gamma - 1)) & | & - ( h W - 1 - \frac{1}{W} (\gamma - 1) \epsilon) & 0 & 1 \end{matrix}\right]$
$\left[\begin{matrix} 1 & \rho W^2 v_j & 0 & | & \frac{1}{W} & 0 & 0 \\ 0 & \rho h W^2 (W^2 v_i v_j + \delta_{ij}) & \gamma \rho W^2 v_i & | & - h W v_i & \delta_{ij} & 0 \\ 0 & \rho W^2 v_j ( h W^2 + (\gamma - 1) \epsilon ) & \rho (\gamma W^2 - (\gamma - 1)) & | & - h W + 1 + \frac{1}{W} (\gamma - 1) \epsilon & 0 & 1 \end{matrix}\right]$
$B \rightarrow \rho (\gamma W^2 - (\gamma - 1)) B - \gamma \rho W^2 v_i C$
$C \rightarrow (\rho (\gamma W^2 - (\gamma - 1)))^{-1} C$
$\left[\begin{matrix} 1 & \rho W^2 v_j & 0 & | & \frac{1}{W} & 0 & 0 \\ 0 & \rho h W^2 (W^2 v_i v_j + \delta_{ij}) \rho (\gamma W^2 - (\gamma - 1)) - \gamma \rho W^2 v_i \rho W^2 v_j ( h W^2 + (\gamma - 1) \epsilon ) & 0 & | & - h W v_i \rho (\gamma W^2 - (\gamma - 1)) - \gamma \rho W^2 v_i (- h W + 1 + \frac{1}{W} (\gamma - 1) \epsilon) & \delta_{ij} \rho (\gamma W^2 - (\gamma - 1)) & - \gamma \rho W^2 v_i \\ 0 & W^2 v_j ( h W^2 + (\gamma - 1) \epsilon ) / (\gamma W^2 - (\gamma - 1)) & 1 & | & (- h W (W - 1) + (\gamma - 1) \epsilon) / (\rho W (\gamma W^2 - (\gamma - 1))) & 0 & (\rho (\gamma W^2 - (\gamma - 1)))^{-1} \end{matrix}\right]$
Let $\Phi = \gamma W^2 - \gamma + 1$
$\left[\begin{matrix} 1 & \rho W^2 v_j & 0 & | & \frac{1}{W} & 0 & 0 \\ 0 & \rho^2 W^2 ( W^2 v_i v_j (\gamma - 1) (1 - 2 h) + h \Phi \delta_{ij} ) & 0 & | & -\rho W v_i (\gamma W - \gamma + 1) & \delta_{ij} \rho \Phi & - \gamma \rho W^2 v_i \\ 0 & W^2 v_j ( W^2 + (\gamma W^2 + \gamma - 1) \epsilon ) / \Phi & 1 & | & ( W (1 - W) h + (\gamma - 1) \epsilon ) / (\rho W \Phi) & 0 & 1 / (\rho \Phi) \end{matrix}\right]$

Finding the inverse of ${A^i}_j = a \delta^i_j + b v^i v_j$
Let ${(A^{-1})^i}_j = c \delta^i_j + d v^i v_j$
${A^i}_k {(A^{-1})^k}_j$
$= (a \delta^i_k + b v^i v_k) (c \delta^k_j + d v^k v_j)$
$= a c \delta^i_j + (a d + b c + b d v^k v_k) v^i v_j$
$= \delta^i_j$ for
$a c = 1$, so $c = \frac{1}{a}$
and
$a d + b c + b d v^k v_k = 0$
$d (a + b v^2) + b c = 0$
$d = -\frac{b}{a (a + b v^2)}$
so
${(A^{-1})^i}_j = \frac{1}{a} \delta^i_j - \frac{b}{a (a + b v^2)} v^i v_j$
$= \frac{1}{a} (\delta^i_j - \frac{b}{a + b v^2} \cdot v^i v_j)$
$= \frac{1}{a} (\delta^i_j - \frac{1}{\frac{a}{b} + v^2} \cdot v^i v_j)$

Let ${\alpha^i}_j = \rho^2 W^2 (h \Phi \delta^i_j + W^2 (\gamma - 1) (1 - 2 h) v^i v_j)$
Therefore ${(\alpha^{-1})^i}_j = \frac{1}{\rho^2 W^2 h \Phi} \cdot ( \delta^i_j - \frac{W^2 (\gamma - 1) (1 - 2 h)}{ h \Phi + W^2 (\gamma - 1) (1 - 2 h) v^2} \cdot v^i v_j )$
$= \frac{1}{\rho^2 W^2 h \Phi} \cdot ( \delta^i_j + \frac{ W^2 (2 h - 1) (\gamma - 1) }{ W^2 + \epsilon \gamma (1 - (\gamma - 2) (W^2 - 1)) } \cdot v^i v_j )$

$ {(\alpha^{-1})^i}_j v^j$
$= \frac{1}{\rho^2 W^2 h \Phi} \cdot ( \delta^i_j - \frac{W^2 (\gamma - 1) (1 - 2 h)}{ h \Phi + W^2 (\gamma - 1) (1 - 2 h) v^2} \cdot v^i v_j ) v^j$
$= \frac{1}{\rho^2 W^2 h \Phi} \cdot ( v^i - \frac{W^2 v^2 (\gamma - 1) (1 - 2 h)}{ h \Phi + W^2 (\gamma - 1) (1 - 2 h) v^2} \cdot v^i )$
$= \frac{1}{\rho^2 W^2 h \Phi} \cdot ( v^i - \frac{(W^2 - 1) (\gamma - 1) (1 - 2 h)}{ h \Phi + (1 - 2 h) (\gamma - 1) (W^2 - 1) } \cdot v^i )$
$= \frac{1}{\rho^2 W^2 h \Phi} \cdot \left( 1 - \frac{1}{1 + \frac{h \Phi}{(W^2 - 1) (\gamma - 1) (1 - 2 h)}} \right) v^i$
aka...
$= \frac{1}{\rho^2 W^2 h \Phi} ( \frac{ h \Phi }{ h \Phi + (1 - 2 h) (\gamma - 1) (W^2 - 1) } ) v^i$
$= \frac{1}{ (\rho^2 W^2)(h \Phi + (1 - 2 h) (\gamma - 1) (W^2 - 1)) } \cdot v^i$

$ {(\alpha^{-1})^i}_j v_i v^j = \frac{1}{ (\rho^2 W^2)(h \Phi + (1 - 2 h) (\gamma - 1) (W^2 - 1)) } \cdot v^2$
$= \frac{W^2 - 1}{ (\rho^2 W^4)(h \Phi + (1 - 2 h) (\gamma - 1) (W^2 - 1)) } $

$B^i \rightarrow {(\alpha^{-1})^i}_k B^k$
$\left[\begin{matrix} 1 & \rho W^2 v_j & 0 & | & \frac{1}{W} & 0 & 0 \\ 0 & \delta^i_j & 0 & | & -\rho W (\gamma W - \gamma + 1) {(\alpha^{-1})^i}_k v^k & \rho \Phi {(\alpha^{-1})^i}_j & -\gamma \rho W^2 {(\alpha^{-1})^i}_k v^k \\ 0 & W^2 v_j ( W^2 + (\gamma W^2 + \gamma - 1) \epsilon ) / \Phi & 1 & | & ( W (1 - W) h + (\gamma - 1) \epsilon ) / (\rho W \Phi) & 0 & 1 / (\rho \Phi) \end{matrix}\right]$
$A^j \rightarrow A^j - {B^i}_j \cdot \rho W^2 v_i$
$C^j \rightarrow C^j - {B^i}_j \cdot W^2 (W^2 + (\gamma W^2 + \gamma - 1) \epsilon) / \Phi \cdot v_i$
$\left[\begin{matrix} 1 & 0 & 0 & | & \frac{1}{W} + \rho W (\gamma W - \gamma + 1) {(\alpha^{-1})^l}_k v^k \cdot \rho W^2 v_l & -\rho \Phi {(\alpha^{-1})^l}_j \cdot \rho W^2 v_l & \gamma \rho W^2 {(\alpha^{-1})^l}_k v^k \cdot \rho W^2 v_l \\ 0 & \delta^i_j & 0 & | & -\rho W (\gamma W - \gamma + 1) {(\alpha^{-1})^i}_k v^k & \rho \Phi {(\alpha^{-1})^i}_j & -\gamma \rho W^2 {(\alpha^{-1})^i}_k v^k \\ 0 & 0 & 1 & | & ( W (1 - W) h + (\gamma - 1) \epsilon ) / (\rho W \Phi) + \rho W (\gamma W - \gamma + 1) {(\alpha^{-1})^l}_k v^k \cdot W^2 (W^2 + (\gamma W^2 + \gamma - 1) \epsilon) / \Phi \cdot v_l & -\rho \Phi {(\alpha^{-1})^l}_j \cdot W^2 (W^2 + (\gamma W^2 + \gamma - 1) \epsilon) / \Phi \cdot v_l & 1 / (\rho \Phi) + \gamma \rho W^2 {(\alpha^{-1})^l}_k v^k \cdot W^2 (W^2 + (\gamma W^2 + \gamma - 1) \epsilon) / \Phi \cdot v_l \end{matrix}\right]$
$\left[\begin{matrix} 1 & 0 & 0 & | & \frac{ 2 h W^2 - \Phi h + \Phi - 1 }{ W (h \Phi + (1 - 2 h) (\Phi - W^2) ) } & \frac{ -\Phi v_j }{ h \Phi + (1 - 2 h) (\Phi - W^2) } & \frac{ \Phi - 1 }{ h \Phi + (1 - 2 h) (\Phi - W^2) } \\ 0 & \delta^i_j & 0 & | & \frac{ -(\gamma W - \gamma + 1) v^i }{ \rho W (h \Phi + (1 - 2 h) (\Phi - W^2)) } & \frac{1}{\rho W^2 h} \cdot ( \delta^i_j + \frac{ W^2 (2 h - 1) (\gamma - 1) }{ W^2 + \epsilon \gamma (1 - (\gamma - 2) (W^2 - 1)) } \cdot v^i v_j ) & \frac{ -\gamma v^i }{ \rho (h \Phi + (1 - 2 h) (\Phi - W^2)) } \\ 0 & 0 & 1 & | & \frac{ h^2 (W (1 - W) + 1) - (\epsilon + 1) h - (\gamma W - \gamma + 1) (W^2 - 1) \epsilon + (h W - h W^2 + h + (W^2 - 1) \epsilon + 1) (1 - 2 h) }{ \rho W (h \Phi + (1 - 2 h) (\Phi - W^2)) } + \frac{ (h W - W^2 + 2) (2 h - 1) W^2 + W^2 (\gamma W - \gamma + 1) (W^2 - 1) (2 h - 1) }{ \Phi \rho W (h \Phi + (1 - 2 h) (\Phi - W^2)) } & \frac{ -(h W^2 + \epsilon (\gamma - 1)) v_j }{ \rho (h \Phi + (1 - 2 h) (\Phi - W^2)) } & \frac{ h W^2 }{ \rho (h \Phi + (1 - 2 h) (\Phi - W^2)) } \end{matrix}\right]$

$\frac{1}{W (h \Phi + (1 - 2 h) (\Phi - W^2))} \cdot \overset{\downarrow I \rightarrow J}{ \left[\begin{matrix} 2 h W^2 - \Phi h + \Phi - 1 & -W \Phi v_j & W (\Phi - 1) \\ - \frac{1}{\rho} (\gamma W - \gamma + 1) v^i & \frac{1}{\rho W^2 h} \cdot W (h \Phi + (1 - 2 h) (\Phi - W^2)) \cdot ( \delta^i_j + \frac{ W^2 (2 h - 1) (\gamma - 1) }{ W^2 + \epsilon \gamma (1 - (\gamma - 2) (W^2 - 1)) } \cdot v^i v_j ) & -\frac{1}{\rho} W \gamma v^i \\ \frac{1}{\Phi \rho} ( + (\gamma - 1) (W^2 - 1) (2 - W) + h ( + \epsilon - 2 W^2 \epsilon - 2 - W^2 + 3 \gamma - W^2 \gamma + W^5 \gamma + W (\gamma - 1) ( - 2 + h - 2 W^3 ) + (\gamma - 2) ( + W^3 + h - h W^2 2 - h W^3 + h W^4 ) ) ) & -\frac{1}{\rho} W (h W^2 + \epsilon (\gamma - 1)) v_j & \frac{1}{\rho} h W^3 \end{matrix}\right] }$

verify
$\frac{\partial U_I}{\partial W_J} \cdot \frac{\partial W_J}{\partial U_K}$

$= \left[\matrix{ W & \rho W^3 v_j & 0 \\ h W^2 v_i & \rho h W^2 (2 W^2 v_i v_j + \delta_{ij}) & \gamma \rho W^2 v_i \\ W (h W - 1) - (\gamma - 1) \epsilon & \rho W^3 v_j (2 h W - 1) & \rho (\gamma W^2 - (\gamma - 1)) }\right] \cdot \frac{1}{W (h \Phi + (1 - 2 h) (\Phi - W^2))} \cdot \left[\begin{matrix} 2 h W^2 - \Phi h + \Phi - 1 & - \frac{1}{\rho} (\gamma W - \gamma + 1) v^k & \frac{1}{\Phi \rho} ( + (\gamma - 1) (W^2 - 1) (2 - W) + h ( + \epsilon - 2 W^2 \epsilon - 2 - W^2 + 3 \gamma - W^2 \gamma + W^5 \gamma + W (\gamma - 1) ( - 2 + h - 2 W^3 ) + (\gamma - 2) ( + W^3 + h - h W^2 2 - h W^3 + h W^4 ) ) ) \\ -W \Phi v_j & \frac{1}{\rho W^2 h} \cdot W (h \Phi + (1 - 2 h) (\Phi - W^2)) \cdot ( \delta^k_j + \frac{ W^2 (2 h - 1) (\gamma - 1) }{ W^2 + \epsilon \gamma (1 - (\gamma - 2) (W^2 - 1)) } \cdot v^k v_j ) & -\frac{1}{\rho} W (h W^2 + \epsilon (\gamma - 1)) v_j \\ W (\Phi - 1) & -\frac{1}{\rho} W \gamma v^k & \frac{1}{\rho} h W^3 \end{matrix}\right]$
$= \frac{1}{W (h \Phi + (1 - 2 h) (\Phi - W^2))} \cdot \overset{\downarrow I \rightarrow K}{\left[\begin{matrix} W (2 h W^2 - \Phi h + \Phi - 1) -W \Phi v_j \rho W^3 v_j & & \\ & & \\ & & \\ \end{matrix}\right]}$
$= \frac{1}{W (h \Phi + (1 - 2 h) (\Phi - W^2))} \cdot \overset{\downarrow I \rightarrow K}{\left[\begin{matrix} W (h \Phi + (1 - 2 h) (\Phi - W^2)) + ... (erroneous terms) ... \gamma W (1 - W) (W^2 - 1) & & & \\ & & \\ & & \\ \end{matrix}\right]}$

change in flux wrt primitive variables:
$\frac{\partial F_{Ij}}{\partial W_K} = \downarrow I \overset{ \rightarrow K }{ \left[\matrix{ D_{,\rho} v_j & D_{,v_k} v_j + D \delta_{jk} & 0 \\ S_{i,\rho} v_j - \delta_{ij} P_{,\rho} & S_{i,v_k} v_j + S_i \delta_{jk} & S_{i,\epsilon} v_j - \delta_{ij} P_{,\epsilon} \\ S_{j,\rho} - D_{,\rho} v_j & S_{j,v_k} - D_{,v_k} v_j - D \delta_{jk} & S_{j,\epsilon} }\right] }$

$\frac{\partial F_{Ij}}{\partial W_K} = \downarrow I \overset{ \rightarrow K }{ \left[\matrix{ W v_j & \rho W^3 v_k v_j + D \delta_{jk} & 0 \\ h W^2 v_i v_j - \delta_{ij} (\gamma - 1) \epsilon & \rho h W^2 (2 W^2 v_i v_k + \delta_{ik}) v_j + S_i \delta_{jk} & \gamma \rho W^2 v_i v_j - \delta_{ij} (\gamma - 1) \rho \\ h W^2 v_j - W v_j & \rho h W^2 (2 W^2 v_j v_k + \delta_{jk}) - \rho W^3 v_k v_j - D \delta_{jk} & \gamma \rho W^2 v_j }\right] }$

$\frac{\partial F_{Ij}}{\partial W_K} = \downarrow I \overset{ \rightarrow K }{ \left[\matrix{ W v_j & D (W^2 v_j v_k + \delta_{jk}) & 0 \\ h W^2 v_i v_j - \delta_{ij} (\gamma - 1) \epsilon & \rho h W^2 (2 W^2 v_i v_k + \delta_{ik}) v_j + S_i \delta_{jk} & \gamma \rho W^2 v_i v_j - \delta_{ij} (\gamma - 1) \rho \\ (h W - 1) W v_j & \rho h W^2 (2 W^2 v_j v_k + \delta_{jk}) - D (W^2 v_j v_k + \delta_{jk}) & \gamma \rho W^2 v_j }\right] }$

acoustic matrix

$\frac{\partial W_I}{\partial U_L} \frac{\partial F_{Lj}}{\partial U_M} \frac{\partial U_M}{\partial W_K}$
...
$= A_{IjK} + \delta_{IK} v_j$

change in flux wrt conservative variables:
$\frac{\partial F_{Ij}}{\partial U_L} = \downarrow I \overset{ \rightarrow K }{ \left[\matrix{ W v_j & D (W^2 v_j v_k + \delta_{jk}) & 0 \\ h W^2 v_i v_j - \delta_{ij} (\gamma - 1) \epsilon & \rho h W^2 (2 W^2 v_i v_k + \delta_{ik}) v_j + S_i \delta_{jk} & \gamma \rho W^2 v_i v_j - \delta_{ij} (\gamma - 1) \rho \\ (h W - 1) W v_j & \rho h W^2 (2 W^2 v_j v_k + \delta_{jk}) - D (W^2 v_j v_k + \delta_{jk}) & \gamma \rho W^2 v_j }\right] } \cdot \downarrow K \overset{ \rightarrow L }{ \left[\matrix{ W & \rho W^3 v_l & 0 \\ h W^2 v_k & \rho h W^2 (2 W^2 v_k v_l + \delta_{kl}) & \gamma \rho W^2 v_k \\ W (h W - 1) - (\gamma - 1) & \rho W^3 v_l (2 h W - 1) & \rho (\gamma W^2 - (\gamma - 1)) }\right]^{-1} }$

Recovering Primitives

Every source, Alcubierre 2008, Marti & Muller 2008, etc, all say the same exact thing.
Except the Mara code uses a paper by Anton & Zanotti http://arxiv.org/pdf/astro-ph/0506063v1.pdf.