The error seems to lie near the end of your derivation, when you claim
$$-\partial_{\mu}(u^{\mu}\rho)=\partial_{t}(v^{0}\rho)-\boldsymbol{\nabla}\cdot\left(\boldsymbol{v}\,\rho\right)$$
(Where I am using bold text to refer to vectors now.) This is, in fact, not the case. Recall that $\partial_{\mu}=(\partial_{t},\boldsymbol{\nabla})$, and requires no minus signs, since derivatives are naturally covariant. On the other hand, $u^{\mu}=(u^{0},\textbf{u})$ naturally transforms like coordinate vectors, and thus naturally has an upstairs index. Thus, the real expression would be
$$-\partial_{\mu}\left(u^{\mu}\rho\right)=-\left(\partial_t(v^{0}\rho)+\boldsymbol{\nabla}\cdot(\boldsymbol{v}\,\rho)\right)$$
With this, the continuity equation is immediately in the form you want it.
The lesson is that, when doing these kinds of calculations, it is important to recall how vectors like $\partial$ and $v$ are defined. Do they naturally have upstairs or downstairs indices? You get negatives when you contract two of the same type of vector, but none when you contract a vector and a covector.
I hope this helps!
Best Answer
Without seeing the manuscript in question, the most obvious reason why is when extremely large ranges of density are expected. If the density varies by orders of magnitude, say in an astrophysics setting, then the log form would keep the numbers in similar scale making it more numerically tractable.
Similar treatment is done for the partial density equations in chemically reacting flows.