Find the general integral of x(y-z) z_x + y(z-x) u_y = z(x-y).
Construct three different solutions z = f(x,y), and describe their domain.

v = {x(y-z), y(z-x), z(x-y)}
u1 = x y z
u2 = x + y + z

gradu1 = D[u1, {{x,y,z}}]
Simplify[gradu1.v]  (* gives 0 iff u1 is a first integral of v *)

Setting x y z -1 gives a solution , z = 1/(xy), which has the domain satisfying x > 0, y > 0.
Note that the graph {(x,y,z) : z = 1/(x y)} is an integral surface of v, and thus a solution surface of the pde.

But u = 1/(x y) - z is NOT a first integral of v.

u = 1/(x y) - z
gradu = D[u, {{x,y,z}}]
graduDotV = Simplify[gradu.v]  (* gives 0 iff u is a first integral of v *)

However, f(x,y) = 1/(x y) does satisfy the PDE because if we evaluate that expression at z = 1/(x y), we get identically 0:

graduDotV/.{z -> 1/(x y)}

The second first integral gives x + y + z = c, or z = f(x,y) = c - x - y, defined on the whole plane.

gradu2 = D[u2, {{x,y,z}}]
Simplify[gradu2.v]  (* gives 0 iff u1 is a first integral of v *)

The set {u1, u2} is functionally independent so we have two FIFIs.

Simplify@Cross[gradu1, gradu2]

For our final solution to the PDE, we can use any F(u1, u2) = 0 and try to solve for z = f(x,y).  If it gets too complicated we can’t solve for z, but here is a family of first integrals.

u = x y z + x + y + z -c
gradu = D[u, {{x,y,z}}]
Simplify[gradu.v]  (* gives 0 iff u1 is a first integral of v *)

Here is an solution surface of the PDE, which is also an integral surface of  V.

ContourPlot3D[x y z +( x +y + z) == 0, {x, -2, 2},{y,-2,2}, {z, -2,2}, AxesLabel->{x,y,z}, ViewPoint->{2.4, -1.3, 2.}]

Note that we can solve for z to get z = f(x,y), which is a solution to the PDE
P(x,y,z) z_x + Q(x,y,z) z_y = R(x,y,z). the domain of the solution is all (x,y) that satisfy 1+ x y > 0
For example x y z + x + y + z = 0 gives z(1+ x y)  + x + y = 0, or z = (- x - y)/(1+ x y).

ContourPlot3D[z == (-x-y)/(1+x y), {x, -2, 2},{y,-2,2}, {z, -2,2}, RegionFunction ->Function[{x,y,z}, 1 + x y >.1],AxesLabel->{x,y,z}, ViewPoint->{2.4, -1.3, 2.}]

Now, the surface z = f(x,y) is an integral surface of the vector field (and a solution surface of the PDE).
I presented to the class the paradox that  u = f(x,y) -z is not an integral surface.

 The resolution to the paradox is that grad u . V = 0 on the surface z = f(x,y).
 In other words, grad u . V is a function of x, y, and z.  It is not identically 0.  But if we substitute in z = f(x,y), we get a function of x and y that is identically 0.
 
Thank you to Shaughn Hirz for figuring this out!
This was a genuine puzzler to me and I wrestled with it for a few days.  As I said in class, the process of resolving paradoxes leads to deeper understanding.  
Try to solve any problem from two different viewpoints. If you get different answers that is the time to wrestle with the paradox.

u = (-x-y)/(1+ x y) -z
gradu = Simplify@D[u, {{x,y,z}}]
graduDotV = Simplify@Expand[gradu.v]  
shouldbe0 = graduDotV/.{z ->(-x-y)/(1+x y)}
Simplify[shouldbe0]
(* gives 0 iff (Simplify works, and u is a first integral of v) *)