Numerically evaluating $\iiint_{C} \frac{ {\rm d}x \, {\rm d}y \, {\rm d}z}{x^2 +y^2 + z^2}$

We can use divergence theorem to say

$$\iiint_{[-1,1]^3} \frac{dV}{r^2} = \iint_{\partial [-1,1]^3} \frac{\hat{r}}{r}\cdot\vec{dS}$$

By symmetry we can choose to only integrate over the square in the $z=1$ plane and evaluate

$$\iint_{[-1,1]^2}\frac{6}{\sqrt{x^2+y^2+1}}\cos\theta\:dA = \iint_{[-1,1]^2}\frac{6}{x^2+y^2+1}\:dA$$

By even more symmetry we can reduce this to $8$ times an integral over a triangle in the first quadrant

$$\int_0^1\int_0^x \frac{48}{x^2+y^2+1}dydx$$

You can show this is equal to

$$\int_0^{\sinh^{-1}(1)}48\tan^{-1}(\tanh t)\:dt$$

by integrating directly or

$$\int_0^{\frac{\pi}{4}}24\log(\sec^2\theta+1)\:d\theta$$

by integrating in polar coordinates. Numerically these integrals both evaluate to $\approx 15.3482$