Divergence theorem, Gauss's theorem $$\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}
$$
