Under which conditions does this PDE have unique solutions Assume that $f:\mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$ is smooth and consider the linear equation $$\mathrm{div}\, (u)(x) = f(x,u(x)),$$ where $u:\mathbb{R}^n \to \mathbb{R}^n $ is a smooth function.
My problem is to find extra conditions so that the equation has locally unique solutions. This conditions should be either new linear equations involving the derivatives of $u$ and boundary conditions.
I know that this is not a very precise question, but I would be very grateful if somebody could provide me a reference in which this kind of systems of PDE are treated, or give me some idea on how to approach this problem.
EXTRA: The origin of this system comes from a problem in differential geometry, where $u$ is an $(n-1)$ form whose exterior derivative is the volume form $f(\cdot,u(\cdot))$. I want to know which extra conditions I need in order to have a well-defined $u$.
 A: The linearization of the operator $Pu(x)=\mathrm{div}(u)(x)-f(x,u(x))$ around $u_{0}$ is the linear operator $(D_{u_{0}}P)v=\mathrm{div}(v)(x)-df(x,u_{0}(x))v$. This linear operator has a surjective symbol, and so the operator $P$ is what is called an underdetermined elliptic operator. If $u_0$ is a function such that at $x_0$ it holds that $Pu_{0}(x_0)=0$, then $u_0$ is called an infinitesimal solution of the equation $Pu=0$ at $x_0$, and there is a local solution for $Pu=0$ around $x_0$ that asymptotically coincides with $u_0$. For reference on this, see for example Chapter 14.3 in Taylor's Partial Differential Equations, Part 3.
Since the system is underdetermined elliptic, in general the local solution will be neither unique nor sufficiently regular. A way to make $P$ into an elliptic operator, hence with certain uniqueness and regularity clauses, is to do what @Willie Wong suggested in the comments and introduce the constraint $u=\nabla w$ for a scalar function $w$. This turns the equation into a second order, nonlinear elliptic equation, and asymptotically the solution ill be unique. In your case, the linearization become $(D_{w_{0}}P)q=\Delta q(x)-df(x,\nabla w_{0}(x))\nabla q(x)$, which by restricting to an appropriate subdomain $\Omega\subset\mathbb{R}^{n}$ and imposing boundary conditions, for instance Dirichlet ones, is invertible for all $w_{0}$. If you restrict $P$ to a certain Banach space, you can then use the implicit function theorem and conclude some invertibility for $P$ when restricted thus, hence uniqueness of a solution to $Pw=0$.
Regularity in this case becomes an issue, but there are regularity theorems for nonlinear elliptic equations which may help. You can read more on this, for example, in Taylor's book.

EDIT: Demanding $u$ to be a gradient is equivalent to demanding the constraint $du=0$ (since we solve locally anyway), where $d$ is the exterior derivative, given in coordinates by $(du)_{ij}=\partial_{i}u_{j}-\partial_{j}u_{i}$. More generally, imposing a condition $du=F_{0}$ for a fixed $F_{0}$ that satisfied $dF_{0}=0$ also turns the equation into an elliptic one: by the Hodge decomposition, every vector field is of the form $u=\nabla w+\delta\omega+\kappa$ where $\kappa$ is a harmonic vector field and $\omega_{ij}$ is a $2$-form such that $(\delta\omega)_{j}=-\sum_{i}\partial_{i}\omega_{ij}$ (we work in a Euclidean setting, so I move from vector fields to differential forms indiscriminately. The arguments are the same in a Riemannian setting by using the musical isomorphism).
Now, demanding $du=F_{0}$ implies that $d\delta\omega=F_{0}$. The condition $dF_{0}=0$ implies that you can solve this equation uniquely for $\omega$, so it effectively insures that $\omega$ is completely determined by the constraint $du=F_{0}$.
As for the other terms in the Hodge decomposition of $u$, $\kappa$ is a finite dimensional obstruction that vanishes if you assume that you solve the equation $Pu=0$ locally. You are left with determining the scalar function $w$. However, since $\mathrm{div}\circ\delta=0$, you return to the situation where you effectively replace $u=\nabla w$, turning the equation $Pu=0$ into an ellitpic second order equation for $w$:
$Pw(x)=\Delta w(x)-f(x,\nabla w(x)+\delta\omega(x))=0$
Note how $\omega$ is fixed, and the only variable is $w$, and so you can deal with this equation in the same manner as my original answer.
It is worth to note that since the Hodge decomposition is valid for a differential form of any order, this disucssion applies also if you replace $u$ with a differential form of any order, with appropriatly replacing $\nabla$ with the exterior derivative $d$ and $\mathrm{div}$ with the dual of $d$, denoted by $\delta$.
A: Here is another way to make the solutions of the equation locally unique:  What you have is a bundle mapping $f$ from the bundle of $(n{-}1)$-forms on $\mathbb{R}^n$ to the bundle of $n$-forms on $\mathbb{R}^n$, and you want to study the equation $\mathrm{d}u = f(u)$, which is underdetermined when $n>1$.
The simplest thing would be to fix a (nonvanishing) $(n{-}1)$-form $\alpha$ and consider $u= v\,\alpha$ where $v$ is a single unknown function.  Then the equation becomes $\mathrm{d}v\wedge\alpha + v\,\mathrm{d}\alpha = f(v\,\alpha)$, which is a single first order PDE for the unknown function $v$.  For example, in local coordinates, if you take $\alpha = \mathrm{d}x^2\wedge\cdots\mathrm{d}x^n$, then the equation has the form
$$
\frac{\partial v}{\partial x^1} = F(x^1,\ldots,x^n,v)
$$
for some function $F$.  Then you can get uniqueness (and local existence) by specifying the values of $v$ along a hypersurface $x^1=const$.
