Alright here's the result. Maybe it should go in my 'Differential Geometry' folder:

$\int_\Omega \nabla \cdot F dV$
$= \int_\Omega \nabla_i F^i (\star 1)$
This is equal to the wedge of all differential forms times the volume of a differential, calculated as the determinant of the change of basis from Cartesian to the curvilinear coordinates:
$= \int_\Omega \nabla_i F^i V dx^1 dx^2 ... dx^n$
The volume of the differential / determinant of the change of basis from Cartesian to curvilinear is also equal to the square root of the determinant of the coordinate metric tensor:
$= \int_\Omega \nabla_i F^i \sqrt{|g|} dx^1 dx^2 ... dx^n$
This is equal to the divergence of the orthonormal component vector using an associated orthonormal non-coordinate (antisymmetric) Levi-Civita covariant derivative:
$= \int_\Omega \nabla_{\hat{i}} F^{\hat{i}} V dx^1 dx^2 \wedge ... \wedge dx^n$
by the Divergence theorem this is equal to:
$= \int_{\partial\Omega} F^{i_1} \epsilon_{i_1 ... i_n} dx^{i_2} \wedge ... \wedge dx^{i_n}$
Where $\epsilon_{i_1 ... i_n}$ is the coordinate-component Levi-Civita permutation tensor, which is weighted by the volume of the differential.
This is equal to the $\pm 1$ permutation (which I will call $\bar{\epsilon}_{i_1 ... i_n}$) times the volume of the differential:
$= \int_{\partial\Omega} F^{i_1} V {\bar{\epsilon}}_{i_1 ... i_n} dx^{i_2} \wedge ... \wedge dx^{i_n}$.
$= \int_{\partial\Omega} F^{i_1} \sqrt{|g|} {\bar{\epsilon}}_{i_1 ... i_n} dx^{i_2} \wedge ... \wedge dx^{i_n}$.
This can be reproduced in an orthonormal non-coordinate system, replacing the coordinate-component Levi-Civita permutation with an orthonormal-basis non-coordinate-component Levi-Civita permutation tensor (weighted by the determinant of the orthonormal non-coordinate metric, which is identity, which has determinant of 1).
$= \int_{\partial\Omega} F^{\hat{i}_1} \epsilon_{\hat{i}_1 ... \hat{i}_n} e^{\hat{i}_2} \wedge ... \wedge e^{\hat{i}_n}$.
$= \int_{\partial\Omega} F^{\hat{i}_1} \sqrt{|\hat{g}|} \bar{\epsilon}_{\hat{i}_1 ... \hat{i}_n} e^{\hat{i}_2} \wedge ... \wedge e^{\hat{i}_n}$.
$= \int_{\partial\Omega} F^{\hat{i}_1} \bar{\epsilon}_{\hat{i}_1 ... \hat{i}_n} e^{\hat{i}_2} \wedge ... \wedge e^{\hat{i}_n}$.
Now we can replace the orthonormal non-coordinate differential form basis with its respective coordinate basis.
$= \int_{\partial\Omega} F^{\hat{i}_1} \bar{\epsilon}_{\hat{i}_1 ... \hat{i}_n} \cdot |e_{i_2}| ... |e_{i_n}| \cdot dx^{i_2} \wedge ... \wedge dx^{i_n}$.
And last, to recover the original coordinate-component formalism, we can replace the vector index non-coordinate component with a coordinate component:
$= \int_{\partial\Omega} |e_{i_1}| F^{i_1} \bar{\epsilon}_{\hat{i}_1 ... \hat{i}_n} \cdot |e_{i_2}| ... |e_{i_n}| \cdot dx^{i_2} \wedge ... \wedge dx^{i_n}$.
$= \int_{\partial\Omega} F^{i_1} \bar{\epsilon}_{\hat{i}_1 ... \hat{i}_n} \cdot |e_{i_1}| ... |e_{i_n}| \cdot dx^{i_2} \wedge ... \wedge dx^{i_n}$.
Now we can assert that $|e_1| \cdot ... \cdot |e_n| = V$. And of course the $\pm 1$ permutation is independent of coordinate basis.
$= \int_{\partial\Omega} F^{i_1} \bar{\epsilon}_{i_1 ... i_n} V dx^{i_2} \wedge ... \wedge dx^{i_n}$.



Application to finite volume:
Recap on Stokes theorem of tensors in index notation:

$\int_\Omega {F^i}_{;i} dV$
$\int_\Omega {F^i}_{;i} V dx^1 dx^2 dx^3$
$= \int_{\partial \Omega} F^i \epsilon_{ijk} dx^j dx^k$ for j,k fixed.
$= \int_{\partial \Omega} F^i \bar{\epsilon}_{ijk} V dx^j dx^k$
$= \int_{\partial \Omega} F^\hat{i} {e^i}_\hat{i} \bar{\epsilon}_{ijk} V dx^j dx^k$

Finite volume update:
$\int_{t \times \Sigma} {T^{\mu\nu}}_{;\nu} dt dx^1 dx^2 dx^3 = 0$
TODO start from curved spacetime.
Until then...
$\int_t \int_\Omega \partial_t (e_{\hat{i}}(x) \textbf{U}^{\hat{i}}(x)) dV dt + \int_t \int_\Omega e_{\hat{i}}(x) \nabla \cdot \textbf{F}^{\hat{i}}(x) dV dt = \int_t \int_\Omega e_{\hat{i}}(x) \textbf{S}^{\hat{i}}(x) dV dt$
$\int_t \int_\Omega \partial_t (e_{\hat{i}}(x_0) {P^\hat{i}}_{\hat{j}}(x,x_0) \textbf{U}^{\hat{j}}(x)) dV dt + \int_t \int_\Omega e_{\hat{i}}(x_0) {P^\hat{i}}_{\hat{j}}(x,x_0) \nabla \cdot \textbf{F}^{\hat{j}}(x) dV dt = \int_t \int_\Omega e_{\hat{i}}(x) \textbf{S}^{\hat{i}}(x) dV dt$
Assume $\partial_t e_{\hat{i}}(x_0) = 0$
$e_{\hat{i}}(x_0) \int_t \int_\Omega \partial_t {P^\hat{i}}_{\hat{j}}(x,x_0) \textbf{U}^{\hat{j}}(x) dV dt + e_{\hat{i}}(x_0) \int_t \int_\Omega {P^\hat{i}}_{\hat{j}}(x,x_0) \nabla \cdot \textbf{F}^{\hat{j}}(x) dV dt = e_{\hat{i}}(x_0) \int_t \int_\Omega {P^{\hat{i}}}_{\hat{j}}(x,x_0) \textbf{S}^{\hat{j}}(x) dV dt$
Let $e_{\hat{i}}(x_0) U^{\hat{i}}(x_0) = e_{\hat{i}}(x_0) \frac{1}{|\Omega|} \int_\Omega {P^{\hat{i}}}_{\hat{j}}(x, x_0) \textbf{U}^{\hat{j}}(x) (dV = V dx^1 dx^2 dx^3)$
Let $e_{\hat{i}}(x_0) S^{\hat{i}}(x_0) = e_{\hat{i}}(x_0) \frac{1}{\Delta t |\Omega|} \int_t \int_\Omega {P^{\hat{i}}}_{\hat{j}}(x, x_0) \textbf{S}^{\hat{j}}(x) dV dt$
$e_{\hat{i}}(x_0) |\Omega| \left( U^{\hat{i}}(t+\Delta t, x_0) - e_{\hat{i}}(x_0) U^{\hat{i}}(t, x_0) \right) + e_{\hat{i}}(x_0) \Delta t \int_\Omega {P^{\hat{i}}}_{\hat{j}}(x, x_0) {\textbf{F}^{\hat{j}\hat{k}}}_{;\hat{k}}(x) dV = \Delta t |\Omega| e_{\hat{i}}(x_0) S^{\hat{i}}(x_0)$
Using the divergence theorem:
$e_{\hat{i}}(x_0) U^{\hat{i}}(t+\Delta t, x_0) = e_{\hat{i}}(x_0) U^{\hat{i}}(t, x_0) + \Delta t \left( - e_{\hat{i}}(x_0) \frac{1}{|\Omega|} \Sigma_I \int_{\partial\Omega_I} {P^{\hat{i}}}_{\hat{j}}(x, x_0) \textbf{F}^{\hat{j}\hat{k}}(x) \epsilon_{{\hat{k}}lm} dx^l dx^m + e_{\hat{i}}(x_0) S^{\hat{i}}(x_0) \right) $ for fixed l,m
$e_{\hat{i}}(x_0) U^{\hat{i}}(t+\Delta t) = e_{\hat{i}}(x_0) U^{\hat{i}}(t) + \Delta t \left( - e_{\hat{i}}(x_0) \frac{1}{|\Omega|} \Sigma_I \left( {P^\hat{i}}_{\hat{m}}(x_I,x_0) \int_{\partial\Omega_I} {P^\hat{m}}_{\hat{j}}(x,x_I) \textbf{F}^{\hat{j}\hat{k}}(x) {e^k}_{\hat{k}}(x) (\epsilon_{klm}(x) = V(x) \bar\epsilon_{klm}) dx^l dx^m \right) + e_{\hat{i}}(x_0) S^{\hat{i}}(x_0) \right) $
Let $e_{\hat{i}}(x_I) F^{\hat{i}}(x_I) = e_{\hat{i}}(x_I) \frac{1}{|\partial\Omega_I|} \int_{\partial\Omega_I} {P^\hat{i}}_{\hat{j}}(x,x_I) \textbf{F}^{\hat{j}\hat{k}}(x) {e^k}_{\hat{k}}(x) V(x) \bar\epsilon_{klm} dx^l dx^m$
Let $n_{\hat{i}} dS = {e^i}_{\hat{i}}(x) V(x) \bar\epsilon_{ijk} dx^j dx^k$ such that $|dS| = |{e^i}_{\hat{i}}(x) V(x)|$ and $n_{\hat{i}}$ is $\pm 1$ orthogonal to $dx^j dx^k$. (the normal should be Cartesian-norm normalized in is lower form?)
Let $e_{\hat{i}}(x_I) F^{\hat{i}}(x_I) = e_{\hat{i}}(x_I) \frac{1}{|\partial\Omega_I|} \int_{\partial\Omega_I} {P^\hat{i}}_{\hat{j}}(x,x_I) \textbf{F}^{\hat{j}\hat{k}}(x) n_{\hat{k}} dS$
$e_{\hat{i}}(x_0) U^{\hat{i}}(t+\Delta t) = e_{\hat{i}}(x_0) U^{\hat{i}}(t) + \Delta t e_{\hat{i}}(x_0) \left( - \frac{1}{|\Omega|} \Sigma_I \left( |\partial\Omega_I| {P^\hat{i}}_{\hat{j}}(x_I,x_0) F^{\hat{j}}(x_I) \right) + S^{\hat{i}}(x_0) \right) $
$U^{\hat{i}}(t+\Delta t) = U^{\hat{i}}(t) + \Delta t \left( - \frac{1}{|\Omega|} \Sigma_I \left( |\partial\Omega_I| {P^\hat{i}}_{\hat{j}}(x_I,x_0) F^{\hat{j}}(x_I) \right) + S^{\hat{i}}(x_0) \right) $

So in spherical coordintes a curvilinear component vector divergence can be exchanged, via Divergence theorem, with the dot of a normal with the vector.
But what about all those extra connection coefficient terms? It seems like our example of $F = e^\hat{\theta}$ absorbed some of these terms: that $\nabla \cdot F$ does include connections but Divergence transforms it into $n \cdot F$ which does not.
So my new challenges:
1) Verify Stokes' theorem one more time for a vector, using something more complex.
2) Verify Stokes' theorem for a (2 0) tensor. Make sure there are no connections.



And here's where it came from:

Divergence Theorem:

$\int_V \nabla \cdot F dV = \int_S n \cdot F dS$

Stokes Theorem:

$\int_V (\nabla \times F) dV = \int_S n \times F dS$

Ok so what is n? What is its magnitude? What norm is used to measure it?



Example:
Spherical coordinates
$\left[\begin{matrix} x \\ y \\ z \end{matrix}\right] = \left[\begin{matrix} r cos(\phi) sin(\theta) \\ r sin(\phi) sin(\theta) \\ r cos(\theta) \end{matrix}\right]$

Inverse:
Let $r_{xy} = \sqrt{x^2 + y^2}$
$r = \sqrt{x^2 + y^2 + z^2} = \sqrt{(r_{xy})^2 + z^2}$
$\theta = tan^{-1}(\frac{r_{xy}}{z})$
$\phi = tan^{-1}(\frac{y}{x})$
$cos(\phi) = \frac{x}{r_{xy}}$
$sin(\phi) = \frac{y}{r_{xy}}$
$cos(\theta) = \frac{z}{r}$
$sin(\theta) = \frac{r_{xy}}{r}$

$\partial_x r = \frac{x}{r}$
$\partial_y r = \frac{y}{r}$
$\partial_z r = \frac{z}{r}$
$\partial_x r_{xy} = \frac{x}{r_{xy}}$
$\partial_y r_{xy} = \frac{y}{r_{xy}}$
$\partial_z r_{xy} = 0$
$e_r = \partial_r = \left[\begin{matrix} e_x & e_y & e_z \end{matrix}\right] \left[\begin{matrix} cos(\phi) sin(\theta) \\ sin(\phi) sin(\theta) \\ cos(\theta) \end{matrix}\right]$

$e_\theta = \partial_\theta = \left[\begin{matrix} e_x & e_y & e_z \end{matrix}\right] \left[\begin{matrix} r cos(\phi) cos(\theta) \\ r sin(\phi) cos(\theta) \\ -r sin(\theta) \end{matrix}\right]$

$e_\phi = \partial_\phi = \left[\begin{matrix} e_x & e_y & e_z \end{matrix}\right] \left[\begin{matrix} -r sin(\theta) sin(\phi) \\ r sin(\theta) cos(\phi) \\ 0 \end{matrix}\right]$

$e_{\hat{r}} = e_r = \left[\begin{matrix} e_x & e_y & e_z \end{matrix}\right] \left[\begin{matrix} cos(\phi) sin(\theta) \\ sin(\phi) sin(\theta) \\ cos(\theta) \end{matrix}\right]$

$e_{\hat{\theta}} = \frac{1}{r} e_\theta = \left[\begin{matrix} e_x & e_y & e_z \end{matrix}\right] \left[\begin{matrix} cos(\phi) cos(\theta) \\ sin(\phi) cos(\theta) \\ -sin(\theta) \end{matrix}\right]$

$e_{\hat{\phi}} = \frac{1}{r sin(\theta)} e_\phi = \left[\begin{matrix} e_x & e_y & e_z \end{matrix}\right] \left[\begin{matrix} -sin(\phi) \\ cos(\phi) \\ 0 \end{matrix}\right]$

$\left[\begin{matrix} e_{\hat{r}} \\ e_{\hat{\theta}} \\ e_{\hat{\phi}} \end{matrix}\right] = \left[\begin{matrix} cos(\phi) sin(\theta) & sin(\phi) sin(\theta) & cos(\theta) \\ cos(\phi) cos(\theta) & sin(\phi) cos(\theta) & -sin(\theta) \\ -sin(\phi) & cos(\phi) & 0 \end{matrix}\right] \left[\begin{matrix} e_x \\ e_y \\ e_z \end{matrix}\right]$

$\left[\begin{matrix} e_x \\ e_y \\ e_z \end{matrix}\right] = \left[\begin{matrix} cos(\phi) sin(\theta) & sin(\phi) sin(\theta) & cos(\theta) \\ cos(\phi) cos(\theta) & sin(\phi) cos(\theta) & -sin(\theta) \\ -sin(\phi) & cos(\phi) & 0 \end{matrix}\right]^{-1} \left[\begin{matrix} e_{\hat{r}} \\ e_{\hat{\theta}} \\ e_{\hat{\phi}} \end{matrix}\right] = \left[\begin{matrix} cos(φ) sin(θ) & cos(φ) cos(θ) & -sin(φ) \\ sin(φ) sin(θ) & sin(φ) cos(θ) & cos(φ) \\ cos(θ) & -sin(θ) & 0 \end{matrix}\right] \left[\begin{matrix} e_{\hat{r}} \\ e_{\hat{\theta}} \\ e_{\hat{\phi}} \end{matrix}\right] $

$e_x = \left[\begin{matrix} e_{\hat{r}} & e_{\hat{\theta}} & e_{\hat{\phi}} \end{matrix}\right] \left[\begin{matrix} cos(φ) sin(θ) \\ cos(φ) cos(θ) \\ -sin(φ) \end{matrix}\right]$

$e_y = \left[\begin{matrix} e_{\hat{r}} & e_{\hat{\theta}} & e_{\hat{\phi}} \end{matrix}\right] \left[\begin{matrix} sin(φ) sin(θ) \\ sin(φ) cos(θ) \\ cos(φ) \end{matrix}\right]$

$e_z = \left[\begin{matrix} e_{\hat{r}} & e_{\hat{\theta}} & e_{\hat{\phi}} \end{matrix}\right] \left[\begin{matrix} cos(θ) \\ -sin(θ) \\ 0 \end{matrix}\right]$

Integrating the divergence across the region $\{ r, \theta, \phi \} \in [\frac{1}{2}, 1] \times [0, \frac{\pi}{2}] \times [0, \frac{\pi}{2}]$.

Where $F = e^{\hat{r}} = \hat{g}^{\hat{r}\hat{r}} e_{\hat{r}} = 1 \cdot \partial_r$

$\int_V \nabla \cdot F dV$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \nabla \cdot F \cdot V d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} (\frac{2}{r} + \partial_r 1) \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} 2 r sin(\theta) d\phi d\theta dr $

$= 2 \int_{r=\frac{1}{2}}^{r=1} r dr \int_{\theta=0}^{\theta=\frac{\pi}{2}} sin(\theta) d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= 2 (\frac{1}{2} r^2 |_{r=\frac{1}{2}}^{r=1} (cos(\theta) |_{\theta=0}^{\theta=\frac{\pi}{2}} 1 |_{\phi=0}^{\phi=\frac{\pi}{2}} $

$= 2 (\frac{1}{2} - \frac{1}{8}) (cos(\frac{\pi}{2}) - cos(0)) \frac{\pi}{2}$

$= -\frac{3 \pi}{8}$

vs

$\int_S n \cdot F dS$

$= \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} e_{\hat{r}} \cdot F r^2 sin(\theta) d\phi d\theta |_{r=1} $
$ + \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -e_{\hat{r}} \cdot F r^2 sin(\theta) d\phi d\theta |_{r=\frac{1}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} e_{\hat{\theta}} \cdot F r^2 sin(\theta) dr d\theta |_{\theta=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -e_{\hat{\theta}} \cdot F r^2 sin(\theta) dr d\theta |_{\theta=0} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=\frac{\pi}{2}}^{\theta=0} e_{\hat{\phi}} \cdot F r^2 sin(\theta) dr d\theta |_{\phi=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=\frac{\pi}{2}}^{\theta=0} -e_{\hat{\phi}} \cdot F r^2 sin(\theta) dr d\theta |_{\phi=0} $

$= \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} e_{\hat{r}} \cdot e^{\hat{r}} r^2 sin(\theta) d\phi d\theta |_{r=1} + \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -e_{\hat{r}} \cdot e^{\hat{r}} r^2 sin(\theta) d\phi d\theta |_{r=\frac{1}{2}} $

$= \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} sin(\theta) d\phi d\theta - \frac{1}{4} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} sin(\theta) d\phi d\theta $

$= \frac{3}{4} \int_{\theta=0}^{\theta=\frac{\pi}{2}} sin(\theta) d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= \frac{3}{4} cos(\theta)|_{\theta=0}^{\theta=\frac{\pi}{2}} 1|_{\phi=0}^{\phi=\frac{\pi}{2}} $

$= \frac{3}{4} \cdot -1 \cdot \frac{\pi}{2}$

$= \frac{-3 \pi}{8}$

CHECK. That worked. That was an easy one, since $e_{\hat{r}} = \partial_r$. Let's try a non-coordinate basis now:

Same but with $F = e^{\hat{\phi}} = \hat{g}^{\hat{\phi}\hat{\phi}} e_{\hat{\phi}} = 1 \cdot \frac{1}{r sin(\theta)} \partial_\phi$.

$\int_V \nabla \cdot F dV$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \nabla \cdot e^{\hat{\phi}} \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \frac{1}{r sin(\theta)} \partial_\phi 1 \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} 0 d\phi d\theta dr $

$= 0$

$\int_S n \cdot F dS$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=\frac{\pi}{2}}^{\theta=0} e_{\hat{\phi}} \cdot e^{\hat{\phi}} r^2 sin(\theta) dr d\theta |_{\phi=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=\frac{\pi}{2}}^{\theta=0} -e_{\hat{\phi}} \cdot e^{\hat{\phi}} r^2 sin(\theta) dr d\theta |_{\phi=0} $

$= (1 - 1)\int_{r=\frac{1}{2}}^{r=1} \int_{\theta=\frac{\pi}{2}}^{\theta=0} e_{\hat{\phi}} \cdot e^{\hat{\phi}} r^2 sin(\theta) dr d\theta $

$= 0$

Still good, but zero is easy to duplicate, especially since the integral is axisymmetric, so there is an equal quantity going in as coming out.

Let's try $F = e^{\hat{\theta}}$. Then you have a net amount in from the $\theta+$ surface and nothing coming out of the $\theta-$ surface (which is a single point at the $z+$ axis.
So $F = e^{\hat{\theta}} = F_{\hat{\theta}} e^{\hat{\theta}}$ for $F_{\hat{\theta}} = F^{\hat{\theta}} = 1$.
And $F = F^\theta e_\theta = F_\theta e^\theta$ for $F^\theta = \frac{1}{r}$ and $F_\theta = r$.

$\int_V \nabla \cdot F dV$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \nabla \cdot e^{\hat{\theta}} \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} (\frac{cos(\theta)}{r sin(\theta)} + \frac{1}{r} \partial_\theta 1) \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} r dr \int_{\theta=0}^{\theta=\frac{\pi}{2}} cos(\theta) d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= \frac{3}{8} \cdot 1 \cdot \frac{\pi}{2} $

$= \frac{3\pi}{16}$

$\int_S n \cdot F dS$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} e_{\hat{\theta}} \cdot e^{\hat{\theta}} r^2 sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -e_{\hat{\theta}} \cdot e^{\hat{\theta}} r^2 sin(\theta) dr d\phi |_{\theta=0} $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^2 sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} $
$ - \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^2 sin(\theta) dr d\phi |_{\theta=0} $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^2 sin(\frac{\pi}{2}) dr d\phi - \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^2 sin(0) dr d\phi $

$= \int_{r=\frac{1}{2}}^{r=1} r^2 dr \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= (\frac{1}{3} r^3 |_{r=\frac{1}{2}}^{r=1} \cdot 1|_{\phi=0}^{\phi=\frac{\pi}{2}} $

$= (\frac{1}{3} - \frac{1}{24}) \frac{\pi}{2}$

$= \frac{7 \pi}{48}$

Alright, finally, something is different.

Let's try performing the same surface integral but with a coordinate normal instead of a non-coordinate normal.

$\int_S n \cdot F dS$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} (e_\theta \cdot e^{\hat{\theta}} = r e_{\hat{\theta}} \cdot e^{\hat{\theta}}) r^2 sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -(e_\theta \cdot e^{\hat{\theta}} = r e_{\hat{\theta}} \cdot e^{\hat{\theta}}) r^2 sin(\theta) dr d\phi |_{\theta=0} $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^3 sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} $
$ - \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^3 sin(\theta) dr d\phi |_{\theta=0} $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r^3 dr d\phi $

$= \int_{r=\frac{1}{2}}^{r=1} r^3 dr \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= (\frac{1}{4} r^4|_{r=\frac{1}{2}}^{r=1} \cdot 1|_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= \frac{15}{64} \cdot \frac{\pi}{2} $

$= \frac{15\pi}{128}$

Still not good. Try again but ...

$\int_S n \cdot F dS$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} (e_\theta \cdot e^{\hat{\theta}} = \frac{1}{r}) r^2 sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -(e_\theta \cdot e^{\hat{\theta}} = \frac{1}{r}) r^2 sin(\theta) dr d\phi |_{\theta=0} $

...but what kind of inner product is this?
(maybe I'm forgetting that $e_{\hat{\theta}} = \frac{1}{r} e_\theta$ means $F^{\hat{\theta}} = r F^\theta$ such that $e_{\hat{\theta}} F^{\hat{\theta}} = e_\theta F^\theta$.

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} (e_\theta \cdot e_{\hat{\theta}} g^{\theta\theta} = r \frac{1}{r^2}?) r^2 sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} $
$ + \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} -(e_\theta \cdot e_{\hat{\theta}} g^{\theta\theta} = r \frac{1}{r^2}?) r^2 sin(\theta) dr d\phi |_{\theta=0} $

Ohh that's right, the normal has to be computed from the Levi-Civita permutation tensor, which is weighted by the (inverse for contravariant) metric determinant.

$\int_S n \cdot F dS$

$= \int_S F^i \epsilon_{ijk} dS^{jk}$

$= \int_{r \times \phi} F^\theta \epsilon_{\theta \phi r} d\phi dr$

$= \int_{r \times \phi} F^\theta \bar{\epsilon}_{\theta \phi r} r^2 sin(\theta) d\phi dr$

$= \int_{r \times \phi} \frac{1}{r} \cdot r^2 sin(\theta) d\phi dr$

$= \int_{r \times \phi} r sin(\theta) d\phi dr$

...and now things look correct.
Mind you, I had to exchange everything into coordinate components. If I was to switch to non-coordinate components, then what would I use for the volume form?
Well, the volume form is $d\phi \wedge dr$, so it should have a magnitude of $(|e_r| = 1) \cdot (|e_\phi| = r sin(\theta)$.
So another way to write this would be:

$= \int_{r \times \phi} F^{\hat{\theta}} \epsilon_{\hat{\theta} \hat{\phi} \hat{r}} e^{\hat{\phi}} \wedge e^{\hat{r}}$

$= \int_{r \times \phi} (F^{\hat{\theta}} = 1) \cdot (\epsilon_{\hat{\theta} \hat{\phi} \hat{r}} = {\bar{\epsilon}}_{\hat{\theta} \hat{\phi} \hat{r}} = 1) \cdot (r sin(\theta) d\phi) \wedge (dr)$

...and once again, things look correct.
On to calculating the correct result:

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r sin(\theta) dr d\phi |_{\theta=\frac{\pi}{2}} - \int_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r sin(\theta) dr d\phi |_{\theta=0} $

$= \int_{r=\frac{1}{2}}^{r=1} r dr \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= (\frac{1}{2} r^2|_{r=\frac{1}{2}}^{r=1} \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= \frac{3}{8} \cdot 1 \cdot \frac{\pi}{2}$

$= \frac{3 \pi}{16}$

I'm going to verify the divergence with Cartesian components, since the equivalent Cartesian component surface integral is not matching this...

$\int_V \nabla \cdot F dV$

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \nabla \cdot e^{\hat{\theta}} \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \left[\begin{matrix} \partial_x & \partial_y & \partial_z \end{matrix}\right] \cdot \left( \left[\begin{matrix} cos(\phi) cos(\theta) \\ sin(\phi) cos(\theta) \\ -sin(\theta) \end{matrix}\right] = \left[\begin{matrix} \frac{x}{r_{xy}} \frac{z}{r} \\ \frac{y}{r_{xy}} \frac{z}{r} \\ -\frac{r_{xy}}{r} \end{matrix}\right] \right) \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \left[\begin{matrix} \partial_x & \partial_y & \partial_z \end{matrix}\right] \cdot \left[\begin{matrix} \frac{x}{r_{xy}} \frac{z}{r} \\ \frac{y}{r_{xy}} \frac{z}{r} \\ -\frac{r_{xy}}{r} \end{matrix}\right] \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} ( \partial_x (\frac{x}{r_{xy}} \frac{z}{r}) + \partial_y (\frac{y}{r_{xy}} \frac{z}{r}) - \partial_z (\frac{r_{xy}}{r}) ) \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} ( z ( \frac{1}{r_{xy}} \frac{1}{r} - x \frac{1}{(r_{xy})^2} \frac{x}{r_{xy}} \frac{1}{r} - x \frac{1}{r_{xy}} \frac{1}{r^2} \frac{x}{r} ) + z ( \frac{1}{r_{xy}} \frac{1}{r} - y \frac{1}{(r_{xy})^2} \frac{y}{r_{xy}} \frac{1}{r} - y \frac{1}{r_{xy}} \frac{1}{r^2} \frac{y}{r} ) + r_{xy} \frac{1}{r^2} \frac{z}{r} ) \cdot r^2 sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} (\frac{z}{r_{xy}} = \frac{cos(\theta)}{sin(\theta)}) \cdot r sin(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} r cos(\theta) d\phi d\theta dr $

$= \int_{r=\frac{1}{2}}^{r=1} r dr \int_{\theta=0}^{\theta=\frac{\pi}{2}} cos(\theta) d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi $

$= \frac{3}{8} \cdot 1 \cdot \frac{\pi}{2} $

$= \frac{3\pi}{16}$

Sure enough, the Cartesian-component divergence integral agrees with the spherical-component divergence integral.



Oops I didn't need this:

$\int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} e_{\hat{r}} \cdot r^2 sin(\theta) d\phi d\theta dr $

$\int_{r=\frac{1}{2}}^{r=1} \int_{\theta=0}^{\theta=\frac{\pi}{2}} \int_{\phi=0}^{\phi=\frac{\pi}{2}} \left[\begin{matrix} cos(\phi) sin(\theta) \\ sin(\phi) sin(\theta) \\ cos(\theta) \end{matrix}\right] \cdot r^2 sin(\theta) d\phi d\theta dr $

$\int_{r=\frac{1}{2}}^{r=1} r^2 dr \cdot \left[\begin{matrix} \int_{\theta=0}^{\theta=\frac{\pi}{2}} sin(\theta)^2 d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} cos(\phi) d\phi \\ \int_{\theta=0}^{\theta=\frac{\pi}{2}} sin(\theta)^2 d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} sin(\phi) d\phi \\ \int_{\theta=0}^{\theta=\frac{\pi}{2}} cos(\theta) sin(\theta) d\theta \int_{\phi=0}^{\phi=\frac{\pi}{2}} d\phi \end{matrix}\right] $

$(1^2 - (\frac{1}{2})^2) \cdot \left[\begin{matrix} ( \frac{1}{2} (\frac{\pi}{2} - sin(\frac{\pi}{2}) cos(\frac{\pi}{2})) - \frac{1}{2} (-sin(0) cos(0)) ) (sin(\frac{\pi}{2}) - sin(0)) \\ ( \frac{1}{2} (\frac{\pi}{2} - sin(\frac{\pi}{2}) cos(\frac{\pi}{2})) - \frac{1}{2} (-sin(0) cos(0)) ) (sin(\frac{\pi}{2}) - sin(0)) \\ ( \frac{1}{2} sin(\frac{\pi}{2})^2 - \frac{1}{2} sin(0)^2 ) \frac{\pi}{2} \end{matrix}\right] $

$\frac{3}{4} \cdot \left[\begin{matrix} \frac{\pi}{4} \\ \frac{\pi}{4} \\ \frac{\pi}{4} \end{matrix}\right] $

$\left[\begin{matrix} \frac{3\pi}{16} \\ \frac{3\pi}{16} \\ \frac{3\pi}{16} \end{matrix}\right] $



Alright, trying for a more complex function, something that can include the connections.
(This means make sure you have components in all directions, especially the r and $\theta$ directions for the sake of connections, and make sure it differentiates in each of the directions.)