Mdthbs

Conservation of charge or rest mass can be written in this way and it is Lorentz invariant
$$
nabla cdot (rho mathbf{u}) + frac{partial rho}{partial t} = 0
$$
So we could be tempted to naively write conservation of energy in this way (I use $gamma_u$ for particle in motion at speed $mathbf{u}$ to not making confusion with $gamma$ relative to speed of $S'$)
$$
nabla cdot (gamma_u rho mathbf{u}) + frac{partial (gamma_u rho)}{partial t} = 0
$$
But this doesn't look Lorentz invariant. I wrong? Exploiting vector identity $nabla cdot (Psi mathbf{A}) = Psi (nabla cdot mathbf{A}) + mathbf{A} cdot (nabla Psi)$ (whit $Psi=gamma_u $) and exploiting conservation of mass, this equation became
$$
left(
mathbf{u} cdot nabla + frac{partial}{partial t} right) gamma_u = 0
$$
where mass is strangely disappeared. But transforming the corresponding primed equation with
$$
frac{partial}{partial x'} = gamma left( frac{partial}{partial x } + frac{v}{c^2} frac{partial}{partial t } right)
$$
$$
frac{partial}{partial y'} = frac{partial}{partial y}
$$
$$
frac{partial}{partial z'} = frac{partial}{partial z}
$$
$$
frac{partial}{partial t'} = gamma left( frac{partial}{partial t } + v frac{partial}{partial x } right)
$$
$$
u_x' = frac{u_x - v}{1-frac{u_x v}{c^2}}
$$
$$
u_y' = frac{u_y}{gamma left( 1-frac{u_x v}{c^2} right)}
$$
$$
u_z' = frac{u_z}{gamma left( 1-frac{u_x v}{c^2} right)}
$$
$$
gamma_{u'} = gamma gamma_u left( 1 - frac{u_x v}{c^2} right)
$$
we get
$$
left(
mathbf{u} cdot nabla + frac{partial}{partial t} right) left[ gamma_u left(1-frac{u_x v}{c^2} right) right] = 0
$$
That it is different than $left(
mathbf{u} cdot nabla + frac{partial}{partial t} right) gamma_u = 0$ written above. Another road could be exploiting
$$
frac{partial gamma_u}{partial x_i} =
frac{gamma_u^3}{c^2} mathbf{u} cdot frac{partial mathbf{u}}{partial x_i} qquad textrm{where $x_i=x,y,z,t$}
$$
to transform $left(
mathbf{u} cdot nabla + frac{partial}{partial t} right) gamma_u = 0$ into
$$
mathbf{u} cdot left(
mathbf{u} cdot nabla + frac{partial}{partial t} right) mathbf{u} = 0
$$
but this equation too doesn't lead to the invariance (although $left(
mathbf{u} cdot nabla + frac{partial}{partial t} right) mathbf{u} = 0$ is actually invariant). There is a way to check the invariance, or writing conservation of energy in that simple way is incorrect?

Your sense of caution is correct. Energy conservation does not work like this, because energy is not a scalar invariant quantity (contrast with electric charge). This means the quantity you wrote as a 4-divergence of energy flux is not a 4-vector. But we can express energy conservation by going one step further into relativity, using the stress-energy tensor. This is a second rank tensor $bf T$ whose components express energy per unit volume, momentum per unit volume, flux of energy, pressure and sheer stress. All of these are involved in a consideration of energy and momentum passing from one place to another, or between one system and another. The conservation of energy and momentum is expressed
$$
partial_mu T^{mu b} = 0
$$
which is a shorthand for
$$
sum_{mu=0}^3 frac{partial}{partial x^mu} T^{mu b} = 0
$$
The physics here is quite involved; this answer is just a small pointer.

You mention in your question something that you refer to as "conservation of mass", but you should note that there is no conservation law for mass, unless you mean the conservation of energy, but then it would be better to call it energy.

In general you have to consider the stress energy tensor. If you want only energy conservation (without the stress and momentum part of the tensor) you can take $partial_nu T^{0 nu}=frac{partial}{partial t} omega + nabla vec{S}/c=0$, where $omega$ is the energy density and $vec{S}$ the energy flow density.

We can't just insert $gamma$ into the equation of continuity, which in this case is a formulation of conservation of rest mass for a free particle:

$$
frac{partial rho}{partial t} + nabla cdot (rho mathbf v) = 0,
$$
and expect the resulting equation will still be valid.

Also, although the above equation has the same form in all inertial frames, this by itself does not imply that the 4-tuple $(rho c, rhomathbf v)$ is a four-vector. In this case it is a 4-vector, similarly to electric current density $j^mu$. But there are other cases where the same kind of equation $partial_mu S^mu = 0$ is valid in all frames, but where $S$ is not a four-vector. Notable example is the Poynting energy density and momentum density 3-vector in matter-free space.

Similar things will happen for matter energy; even if (and that is a big if) one could derive such simple equation for this energy, this wouldn't imply the energy 4-tuple is a 4-vector. In fact, in EM theory based on Maxwell's equations there is no way to formulate energy conservation where energy density of matter or EM field is a part of some 4-vector field; one must one 4-tensors of 2nd rank (which are represented by 4x4 entries).