Suppose we wish to find a function $u(\vec x,t)$ on some domain $\Omega\times[0,\infty)$ satisfying the PDE $$u_t=\|\nabla u\|_2^2,$$ where $\Omega\subset\mathbb R^n$, and the gradient $\nabla$ is in the $\vec x$ independent variable. We are given the initial information $u(\vec x,0)=u_0(\vec x)$ for a prescribed function $u_0:\Omega\to\mathbb R$, as well as Neumann boundary conditions $\nabla u(\vec x,t)\cdot\vec n(\vec x)=0$ for all $\vec x\in\partial \Omega$.
Can we find a closed-form solution to this PDE, e.g. using the method of characteristics?
Following the second example in this post as well as these course notes, it seems like the characteristic curves of this PDE satisfy: $$ \begin{array}{rl} u(t) &=-\|\vec c\|_2^2 t + D\\ \vec x(t) &=-2\vec c t + E, \end{array} $$ but I'm unclear on how to turn these into a solution to the PDE, and whether it's possible to incorporate the Neumann conditions. One can read off a solution $u(\vec x,t)=\|\vec c\|_2^2 t + \vec x \cdot\vec c + D$ for a constant $\vec c$, but this does not satisfy arbitrary initial conditions or the Neumann conditions.
I did find this post, which gives an expression that is close to a closed-form solution via the Hopf-Lax formula for a similar (but not identical) PDE. Perhaps this is as close as we can get?