$$\mathbf F(x,y,z)=x\mathbf i + y\mathbf j+z\mathbf k$$ and $$ D=\left\{x,y,z\in \mathbf R^3 : 0\le z \le1-x^2-y^2\right\}.$$
I want to calculate both
$$\iiint_D\nabla\cdot\mathbf F\,dV = \oint_{\delta D}\mathbf F \cdot\mathbf N\,dS.$$
So far:
$$\nabla f=1+1+1=3$$
$$\oint_{\delta D}f\cdot n\,dS=\iiint_D=3\,dV$$
and after that
$$v(D)=\int_0^\pi\int_0^{2\pi}\int_0^rr^2\sin(\phi)\,dr\,d\theta \,d\phi=\frac{4\pi r^3}{3}.$$
But what is the $r$ in there? Is it $r=1$?
$$\frac{4\pi1^3}{3}=\frac{4\pi}{3}?$$
And after that:
$$\oint_D\mathbf F\cdot \mathbf N\,dS=\int_0^\pi\int_0^{2\pi}r(\theta,\phi)\cdot \mathbf n\,dS=\int_0^\pi\int_0^{2\pi}r^3\sin(\phi)\,d\theta\,d\phi=4\pi r^3,$$ and this $4\pi r^3$ should be equal to $\dfrac{4\pi}{3}$ but is it not. There must be some mistake somewhere but where?
As I said in the comments, $D$ is not a ball. It's bounded below by the unit disk in the $xy$-plane, and above by the paraboloid $z=1-x^2-y^2$. (You may be thinking of $z = \sqrt{1-x^2-y^2}$; that surface is the upper unit hemisphere.) In cylindrical coordinates, $$ D = \left\{(r,\theta,z) : 0 \leq r \leq 1,\ 0 \leq z \leq 1-r^2 \right\} $$ So \begin{align*} \iiint\limits_{D} \nabla \cdot \mathbf{F}\,dV &= \int_0^{2\pi} \int_0^1 \int_0^{1-r^2} 3 r\,dz\,dr \,d\theta \\&= 3\cdot 2\pi \int_0^1 r(1-r^2) \,dr \\&= 3 \cdot 2\pi \cdot \frac{1}{4} = \frac{3\pi}{2} \end{align*}
The boundary of $D$ is the union of two surfaces: \begin{align*} S_1 &= \left\{(r,\theta,0) : 0 \leq r \leq 1 \right\} \\ S_2 &= \left\{(r,\theta,1-r^2) : 0 \leq r \leq 1 \right\} \\ \end{align*} As oriented surfaces, $\partial D = S_2 - S_1$. On $S_1$, the normal vector is vertical, and $\mathbf{F}$ has no vertical component. Therefore, $$ \iint\limits_{S_1} \mathbf{F}\cdot \mathbf{N}\,dS = \iint\limits_{S_1} (x \mathbf{i} + y \mathbf{j}+z\mathbf{k}) \cdot \mathbf{k}\,dA = \iint\limits_{S_1} 0\,dA = 0 $$ On top, $S_2$ is the graph of $z=1-x^2-y^2$, so the vector $\left<2x,2y,1\right>$ is normal to $S_2$. Therefore \begin{align*} \iint\limits_{S_2} \mathbf{F}\cdot \mathbf{N}\,dS &= \iint\limits_{S_1}\left<x,y,1-x^2-y^2\right>\cdot\left<2x,2y,1\right>\,dA \\&= \iint\limits_{S_1}(1+x^2+y^2)\,dA = \iint\limits_{S_1}(1+r^2)\,dA \\&= \int_0^{2\pi} \int_0^1 (1+r^2)r\,dr\,d\theta \\&= 2\pi \int_0^1 \left(r + r^3\right)\,dr = 2\pi \cdot\frac{3}{4} = \frac{3\pi}{2} \end{align*} Therefore $$ \iint\limits_{\partial D} \mathbf{F}\cdot \mathbf{N}\,dS = \frac{3\pi}{2} $$