Calculating Divergence for a Field in Spherical Coordinates Calculate the flux of the following:
F $=\sqrt{x^2+y^2+z^2}⟨ x,y,z⟩$ 
$D = $ the region $1 ≤ x^2 + y^2 + z^2 ≤ 2$ 
My Work So Far
This is obviously supposed to use the Divergence Theorem
$Flux = \int\int_S$F$\cdot$n$dσ$ $= \int\int\int_D div$ F $dV$
With the boundaries of the search in spherical coordinates: $ρ_{1 \to \sqrt2}, \theta_{0\to 2\pi}, \phi_{0 \to \pi}$
The problem now comes from trying to determine $div$ F, which leads to extremely tedious derivatives (doable, but clearly not intended). If I convert F to spherical coordinates immediately, though, it becomes much cleaner:
F $=\rho⟨\rho sin\phi cos\theta,\rho sin\phi sin\theta,\rho cos\phi⟩$ $\to$
F $=⟨\rho^2 sin\phi cos\theta,\rho^2 sin\phi sin\theta,\rho^2 cos\phi⟩$
Great, much better. The problem is, I now don't see a way to calculate the divergence. Because it takes the form:
$div$ F $= \frac{\partial M}{\partial x} + \frac{\partial N}{\partial y} + \frac{\partial P}{\partial z}$ ($M$ being $\rho^2 sin\phi cos\theta$, etc), and there's no longer an $x,y,z$ to take the partial with respect to, it doesn't seem possible to calculate the divergence in spherical coordinates. But it seems like there absolutely should be, and the question is implying that it should be done that way. What am I overlooking, and how could I calculate $div$ F after doing the transformation above?
 A: The divergence in spherical coordinates is given by 
$$\nabla \cdot \vec F=\frac{1}{r^2}\frac{\partial (r^2F_r)}{\partial r}+\frac{1}{r\sin(\theta)}\frac{\partial (\sin(\theta)F_\theta)}{\partial\theta}+\frac{1}{r\sin(\theta)}\frac{\partial F_\phi}{\partial\theta} \tag 1$$
Here, we have $\vec F=\hat r r^2$. Therefore, 
$$\nabla \cdot \vec F=\frac{1}{r^2}\frac{\partial (r^4)}{\partial r}=4r\tag 2$$

NOTE:
If we are unaware of $(1)$, we can use the product rule $\nabla \cdot (\phi \vec A)=\phi\nabla \cdot \vec A+\vec A\cdot \nabla \phi$ with $\phi=r$ and $\vec A=\vec r$.  Then, $\nabla \cdot (\phi \vec A)=4r$ as expected!

Finally, using $(2)$ reveals
$$\begin{align}
\oint_S \vec F\cdot \hat n\,dS&=\int_V \nabla \cdot \vec F\,dV\\\\
&=\int_0^{2\pi}\int_0^\pi \int_1^{\sqrt 2}(4r)\,r^2\sin(\theta)\,dr\,d\theta\,d\phi\\\\
&=12\pi
\end{align}$$

Note that we could have proceeded by directly evaluating the closed surface integral.  To do this, we simply write
$$\begin{align}
\oint_S \vec F\cdot \hat n\,dS&=\int_0^{2\pi}\int_0^\pi \left.\left(r^4\hat r\right)\right|_{r=\sqrt 2}\cdot \hat r\,\sin(\theta)\,d\theta\,d\phi+\int_0^{2\pi}\int_0^\pi \left.\left(r^4\hat r\right)\right|_{r=1}\cdot (-\hat r)\,\sin(\theta)\,d\theta\,d\phi\\\\
&=12 \pi
\end{align}$$
