Is it possible to reverse a gradient ($\vec{\nabla}$) operation?

Yes, it is possible. I will give an example to demonstrate the general procedure. Consider $f(x, y) = x^3 - 3xy^2 + x^2 + y^2 + \log x$. Then $\nabla f = \dfrac{\partial f}{\partial x}\hat{i} + \dfrac{\partial f}{\partial y}\hat{j} = \left(3x^2 - 3y^2 + 2x + \dfrac{1}{x}\right)\hat{i} + (2y -6xy)\hat{j}$.

To reverse this, we look at each component individually. We know that
$\dfrac{\partial f}{\partial x} = 3x^2 - 3y^2 + 2x + \dfrac{1}{x}\\ \dfrac{\partial f}{\partial y} = 2y - 6xy $

Therefore:

$\begin{align} \displaystyle f(x,y) & = \int \dfrac{\partial f}{\partial x}\, dx\\ & = \int 3x^2 - 3y^2 + 2x + \dfrac{1}{x}\, dx\\ & = x^3 - 3xy^2 + x^2 + \log x + u(y) \end{align}$

What is that $u(y)$? It's the "constant" of integration, of course.When we differentiate $f$ with respect to $x$ partially, any term of $f$ not containing $x$ is a constant - this includes terms containing only $y$.

Now to determine $u(y)$, we look at $\dfrac{\partial f}{\partial y}$. We could integrate this with respect to $y$, similar to what we did with $\dfrac{\partial f}{\partial x}$. Then some of the terms after the integration will be common. The terms that are not common are those that constitute $u(y)$. So we need to look for terms of $\dfrac{\partial f}{\partial y}$ that do not contain $x$, and integrate them. Here, $\dfrac{\partial f}{\partial y} = 2y - 6xy$, and the only term not containing $x$ is $2y$. Therefore:

$\displaystyle u(y) = \int 2y\, dy = y^2 + C$.

Thus, $\boxed{f(x, y) = x^3 - 3xy^2 + x^2 + \log x + y^2 + C}$.

In general: $$f(x, y) = \int \dfrac{\partial f}{\partial x} dx + \int \left[\text{terms of $\dfrac{\partial f}{\partial y}$ that do not contain $x$}\right]\, dy$$

I prefer a simpler one-step solution: Choose the curve $\vec\gamma(t) = t\vec x$ on $t\in[0,1]$, with endpoints $\vec\gamma(0)=\vec 0$ and $\vec\gamma(1) = \vec x$, and then the gradient theorem implies that $$\begin{align} f(\vec x) &= f(\vec 0) + \int_0^1 \nabla f(\vec\gamma(t))\cdot\vec\gamma'(t)\,\mathrm dt \\ &= f(\vec 0) + \int_0^1 \nabla f(t\vec x)\cdot\vec x\,\mathrm dt. \end{align}$$ The unknown $f(\vec 0)$ plays the role of the constant of integration.

For example, if we have $\nabla f(x,y)=\begin{bmatrix}3x^2 - 3y^2 + 2x\\2y - 6xy\end{bmatrix}$, then $$\nabla f(tx,ty)\cdot(x,y) = 3 t^2 (x^3-3xy^2)+2 t (x^2+y^2)$$ from which we easily get $$f(x,y)=f(0,0)+x^3-3xy+x^2+y^2.$$ The $\log x$ term in M. Vinay's answer complicates things a little as the function is undefined at $\vec 0$. It is straightforward to get around this by taking a different starting point, say $\vec x_0=(1,0)$ and choosing the curve $\vec\gamma(t)=\vec x_0+t(\vec x-\vec x_0)$.

M. Vinay's answer is good. For more than 2 variables, the process can be chained. For $$F(x,y,z) = f_x x + f_y y + f_z z = \nabla f$$ we can integrate with respect to $x$

$$f(x,y,z)=\int f_x dx + g(y, z)$$

where $g$ is a function only in terms of $y$ and $z$. Then differentiating with respect to $y$

$$f_y(x,y,z)=\frac{\partial}{\partial y} \int f_x dx \ + g_y(y,z) $$

The first term is common and we are left with $g_y(y,z)$, which can be integrated to get $$g(y,z) = \int g_y(y,z) dy + h(z)$$

$$f(x,y,z)=\int f_x dx + \int g_y(y,z) dy + h(z)$$

Then repeating the process by differentiating $f$ with respect to $z$ and then integrating, we get all the terms of $f$:

$$f(x,y,z)=\int f_x (x,y,z) dx + \int g_y(y,z) dy + \int h_z(z) dz$$