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.