You could ask the same question about a region in the plane: can you integrate the divergence of a vector field on a region to obtain a flux through the boundary curve? The answer is yes. If $R$ is any region in the plane with counterclockwise boundary curve $\partial R$, then
$$
\int_{R}\left(\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial y}\right) dA \;=\; \oint_{\partial R} \left(-Q\,dx + P\, dy\right).
$$
for any vector field $P\hspace{0.5pt}\textbf{i}+Q\textbf{j}$. However, this two-dimensional version of the Divergence Theorem isn't anything new: it's just Green's Theorem applied to the vector field $(-Q,P)$. Note that $(-Q,P)$ is obtained from $(P,Q)$ by rotating each vector counterclockwise by $90^\circ$, so the two-dimesnonal Divergence Theorem is just a "rotated" version of Green's Theorem.
The Divergence Theorem for a surface in $\mathbb{R}^3$ is similar, except that it's basically a "rotated" version of the Curl Theorem (sometimes called Stokes' Theorem). Specifically, if $S$ is a surface with unit normal $\textbf{n}$ and $\textbf{F}$ is a vector field tangent to $S$, then we can rotate $\textbf{F}$ counterclockwise $90^\circ$ by taking its cross product with $\textbf{n}$ (so the rotated field is $\textbf{n}\times\textbf{F}$). Then we can define the "divergence" of $\textbf{F}$ on $S$ by
$$
\mathrm{div}_S(\textbf{F}) \;=\; \textbf{n}\cdot \mathrm{curl}(\textbf{n}\times\textbf{F}).
$$
This formula makes sense even if $\textbf{F}$ isn't tangent to $S$, since it ignores any component of $\textbf{F}$ in the normal direction.
The curl theorem tells us that
$$
\int_S \textbf{n}\cdot\mathrm{curl}(\textbf{G})\;dA \;=\; \oint_{\partial S} \textbf{t}\cdot \textbf{G}\; ds
$$
for any vector field $\textbf{G}$, where $\textbf{t}$ is the unit tangent vector to the curve $\partial S$ (oriented using the right-hand rule). Substituting $\textbf{G} = \textbf{n}\times\textbf{F}$ gives
$$
\int_S \mathrm{div}_S(\textbf{F})\,dA \;=\; \oint_{\partial S} \textbf{t}\cdot(\textbf{n}\times\textbf{F})\; ds.
$$
This is the Divergence Theorem on a surface that you're looking for. The triple product $\textbf{t}\cdot(\textbf{n}\times\textbf{F})$ computes the flux of $\textbf{F}$ through the boundary curve. Perhaps a better way to write the same formula is
$$
\int_S \mathrm{div}_S(\textbf{F})\,dA \;=\; \oint_{\partial S} (\textbf{t}\times\textbf{n})\cdot\textbf{F}\; ds,
$$
where $\textbf{t}\times\textbf{n}$ is a unit vector tangent to the surface, perpendicular to the boundary curve, and pointing directly "out".
First off, you mean to say $f(r) = \hat r / r^2$ since as you've noted, you know the divergence is zero when $r \neq 0$ since $$\nabla \cdot f(r) \hat r = \frac{1}{r^2} \frac{\partial }{\partial r} ( r^2 f(r) )$$
but let's try to under stand what happens when $r \to 0$. Since the function doesn't make sense here, let's try to understand this in a weak sense (or the sense of distributions, a surface integral will come about only when we start trying to make sense of this). Personally, I like to notice that $ -\nabla \frac{1}{r} = f(r)$, i.e. $\Delta \frac{1}{r}$ is the same as $\nabla \cdot f(r)$ (where $\Delta$ is the Laplacian operator). Let $\phi \in C^\infty$ be a nice test function, and consider the integral over the ball $B = \{ x \in \mathbb{R}^3 : ||x|| < \epsilon \}$.
$$ \int _B \phi(x) \Delta \left ( \frac{1}{r} \right ) dx $$
Using Green's identity ( i.e . integration by parts), you may check that
$$ \int _B \left [ \phi(x) \Delta \left ( \frac{1}{r} \right ) - \frac{1}{r} \Delta \phi \right ]dx = \oint_{\partial B} \left [ \phi \frac{\partial }{\partial r} \frac{1}{r} - \frac{1}{r} \frac{\partial \phi }{ \partial r} \right ] dS$$
Note that $dx = r^2 \sin \theta dr \, d\theta \, d \varphi$, thus
$$ \int _B \phi(x) \Delta \left ( \frac{1}{r} \right ) dx= \underbrace{ \int_B \frac{1}{r} \Delta \phi dx }_{ \mathcal{O}(\epsilon)}+ \oint_{\partial B} \left [ \underbrace{\phi \frac{\partial }{\partial r} \frac{1}{r} }_{\mathcal{O}(1)}- \underbrace{\frac{1}{r} \frac{\partial \phi }{ \partial r}}_{\mathcal{O}(\epsilon)} \right ] dS$$
Now if we take the limit as the ball shrinks to the origin (the bad point) $\epsilon \to 0$, we see terms of order $\epsilon$ die, and we're left with
$$ \lim_{\epsilon \to 0} \int _B \phi(x) \Delta \left ( \frac{1}{r} \right ) dx = \lim_{ \epsilon \to 0} \oint_{\partial B} \phi \frac{\partial }{\partial r} \frac{1}{r} dS = \lim_{ \epsilon \to 0} -\int_0^\pi \int_0^{2 \pi } \phi(x) \sin \theta d \theta d \varphi = - 4 \pi \phi(0)$$
Thus some people like to write
$$ \Delta \frac{1}{r} = - 4 \pi \delta (r) $$
Best Answer
The answer is quite simple. If I understand correctly, then you are solving the problem of a magnetic field around a infinite conductor with current. Writing down the equation for the magnetic field circulation, you got the $\phi$-th component of the magnetic field (only it is nonzero), i.e. \begin{equation} B_{\phi} = \frac{\mu_0 I}{2\pi r},\\ B_{z} = 0,\\ B_{r} = 0. \end{equation} The divergence in cylindrical coordinates is as follows \begin{equation} \vec{\nabla}\cdot\vec{B} =\frac{1}{r}\frac{\partial(r B_{r})}{\partial r} + \frac{\partial B_z}{\partial z} + \frac{1}{r}\frac{\partial B_{\phi}}{\partial\phi}, \end{equation} which is in your case is \begin{equation} \vec{\nabla}\cdot\left(\frac{\mu_0 I}{2\pi r}\vec{e}_{\phi}\right) = \frac{1}{r}\frac{\partial\left(\frac{\mu_0 I}{2\pi r}\right)}{\partial\phi} = 0, \end{equation} where $\vec{e}_{\phi}$ is the unit vector along $\phi$.