Although I have detailed a solution in the comments, I think it is probably good to provide the full solution here with an explanation.
Eliminating $dt$ in your system of ODEs yields the following form
$$\frac{dx}{x(y-u)} = \frac{dy}{y(u-x)} = \frac{du}{u(x-y)}$$
Noting that $dA/A = d \ln A$, this can be simplified to
$$\frac{d \ln x}{y-u} = \frac{d \ln y}{u-x} = \frac{d \ln u}{x-y}$$
Now, what the expression above is saying is that the vectors $(d \ln x, d \ln y, d \ln u)$ and $(y - u, u - x, x - y)$ are proportional to eachother i.e
$$(d \ln x, d \ln y, d \ln u) \propto (y - u, u - x, x - y)$$
You might notice that the RHS can be written as the curl of two vectors
$$(y - u, u - x, x - y) = (x, y, u) \times (1, 1, 1)$$
so what we really have is
$$(d \ln x, d \ln y, d \ln u) \propto (x, y, u) \times (1, 1, 1)$$
Remembering that the cross product of two vectors $a$ and $b$ yields a third vector $c$ which is orthogonal to both $a$ and $b$, the expression above really says that the tangent of the log of the solution curve is orthogonal to the vectors $(x, y, u)$ and $(1, 1, 1)$. So, if we take the dot product of the LHS with both vectors $(x, y, u)$ and $(1, 1, 1)$, we can get our integral curves
\begin{align}
(d \ln x, d \ln y, d \ln u) \cdot (x, y, u) &= ((x, y, u) \times (1, 1, 1)) \cdot (x, y, u) \\
&= 0 \\
\implies x d \ln x + y d \ln y + u d \ln u &= dx + dy + du \\
&= 0 \\
\implies x + y + u &= c_{1} \\\\
(d \ln x, d \ln y, d \ln u) \cdot (1, 1, 1) &= ((x, y, u) \times (1, 1, 1)) \cdot (1, 1, 1) \\
&= 0 \\
\implies d \ln x + d \ln y + d \ln u &= 0 \\
\implies \ln x y u &= c_{2} \\\\
\therefore \ln x y u &= f(c_{1}) \\
&= f(x + y + u)
\end{align}
and from here the result follows by applying the initial conditions as stated in the comments.
Note that this approach works well here because of the symmetry in the coefficients of the original PDE. Often times it can be quite difficult to determine the form of the two vectors whose curl is proportional to the tangent vector. However, I think this method gives a nice geometric view of what is happening.
Consider fixed $c_1$ and $c_2$, then these two equations determine a curve. Now, let $c_1$ run free but remaining $c_2$ fixed, obviously we have a surface and it is a solution of the PDE. This gives us a tip about how the particular solutions emerge (we got one!). We can do better, we can move $c_1$ and simultaneously move $c_2$, we draw a surface that is solution too. But the simultaneous movement of $c_1$ and $c_2$ is what we call "function" and because this function is not still determined, we say that it is arbitrary. So $c_2=f(c_1)$ or $\psi=f(\phi)$
A particular solution determines $f$. Consider this example: we know the value of $u$ along the line $y=0$, for each $x$ we know $u$, i. e. $u$ is a perfectly known function $g$ of $x$ $u=g(x)$ along $y=0$ (you can imagine $x^2$ or $e^x$). Then each of the curves have to satisfy that requirement, so is, $\phi(x,0,g(x))=c_1$ and $\psi(x,0.g(x))=c_2$. Now, $x$ can be eliminated determining the needed functional relation between $c_1$ and $c_2$
Best Answer
As per usual set-up, we are looking at $\mathbb{R}^3$ coordinates $(x,y,z)$ for surfaces $z(x,y)=u(x,y)$. Geometrically the PDE says the vector field $\mathbf{v}(x,y,z)=(a(x,y,u(x,y)),b(x,y,u(x,y)),c(x,y,u(x,y)))$ is orthogonal to $(u_x(x,y,u(x,y)),u_y(x,y,u(x,y)),-1)$.
The level sets of $\phi$ (and similarly $\psi$) are therefore integral surfaces of the vector field $\mathbf{v}$ in $\mathbb{R}^3$ from $\frac{\mathrm{d}x}{a}=\frac{\mathrm{d}y}{b}=\frac{\mathrm{d}u}{c}$. We want $\phi$ and $\psi$ to have the property that $\phi=C_1$ intersect $\psi=C_2$ transversely for $C_1,C_2\in\mathbb{R}$.
Now $f(v,w)=0$ is (assuming regularity conditions) a smooth curve in $vw$-plane. So $f(\phi,\psi)=0$ defines a smooth union of characteristic curves, so is also an integral surface of $\mathbf{v}$.