Modifying the problem. Rather than consider the PDE $u_{xt}+u\ u_{xx}+\frac{1}{2}u_{x}^2=0$ with initial condition $u(x,0)=u_0(x)$ as asked above, I will consider the following variant.
$$
\text{Solve }u_{xy}+u\ u_{xx} + u_x^2=0\text{ subject to }u(x,0)=f(x).\qquad(\star)
$$
There are three differences between this question and that which was asked originally.
- The coefficient of $u_x^2$ has changed from $\frac{1}{2}$ to $1$.
- The variable $t$ has been renamed to $y$.
- The initial function $u_0(x)$ has been renamed to $f(x)$.
Only (1) represents a significant modification of the problem. It makes the solution more tractable and enables it to be found using an elementary application of the method of characteristics. For these reasons, it is conceivable that this was the intended question.
Note: I will not delve into regularity of the solutions in this answer.
Reduction to a first order quasilinear PDE. Write the equation as
$$
\frac{\partial}{\partial x}\left(u_y+u\ u_x\right)=0.
$$
Thus $(\star)$ is equivalent to
$$
u_y+u\ u_x=g(y),\qquad u(x,0)=f(x),\qquad (\star\star)
$$
where $g(y)$ is an arbitrary function of $y$ (with sufficient regularity).
Method of characteristics.
Perhaps the simplest formulation of the method of characteristics is for quasilinear first order PDEs. These are PDEs of the form
$$a(x,y,u)u_x+b(x,y,u)u_y=c(x,y,u).$$
To solve this equation, one regards the solution as a surface $z=u(x,y)$ in $xyz$-space. Let $s$ parametrize the initial curve $\bigl(s,0,f(s)\bigr)$ and let $t$ be a second parameter, which can be thought of as the distance flowed along a characteristic curve emanating from $\bigl(s,0,f(s)\bigr)$.
The characteristic equations are then
$$
\frac{dx}{dt}=a(x,y,z),\quad \frac{dy}{dt}=b(x,y,z),\quad \frac{dz}{dt}=c(x,y,z).
$$
Returning to our equation $(\star\star)$, this reduces to $a(x,y,u)=u$ and $b(x,y,u)=1$ and $c(x,y,u)=g(y)$.
Thus
$$
\frac{dx}{dt}=z,\quad \frac{dy}{dt}=1,\quad \frac{dz}{dt}=g(y)
$$
with initial conditions $x(0)=s$ and $y(0)=0$ and $z(0)=f(s)$.
The solution to this system is
$$
x=s+zt,\quad y=t,\quad z=f(s)+h(t),
$$
where $h(t)$ is the antiderivative of $g(t)$ satisfying $h(0)=0$. Since $g$ was arbitrary, so is $h$ given $h(0)=0$.
The solution. Now we eliminate all occurrences of $t$ by replacing them with $y$, then eliminate $s$ by writing $s=x-zy$. Finally, replace $z$ with $u$ to obtain the implicit equation
$$
\boxed{u=f(x-uy)+h(y)},
$$
where $h(y)$ is any sufficiently regular function satisfying $h(0)=0$. This is an implicit equation for the general solution of $(\star)$.
TL;DR. Change the $\frac{1}{2}$ in the original question to $1$ to obtain a PDE solvable by the method of characteristics.
The characteristics are the curves $x=x_0+\alpha u|u|^{\alpha-2}t$ along which $u=u(x_0,0)=F(x_0)$ is constant.
With the given initial data, we have
$$
u(x,t) = \left\lbrace
\begin{aligned}
&0 &&\text{if}\quad x <0\\
&U (x/t) &&\text{if}\quad 0 \leq x \leq \alpha t\\
&1 &&\text{if}\quad \alpha t <x
\end{aligned}
\right.
$$
where $U (x/t)$ is a rarefaction wave solution. The latter is obtained from assuming a smooth solution of the form $u (x,t) = U (\xi)$ with $\xi=x/t$. Using this self-similarity Ansatz and the PDE, one obtains
$$
U (x/t) = \left.\left (\xi\to\alpha \xi|\xi|^{\alpha-2}\right)^{-1}\right|_{\xi=x/t} = \left(\frac{x}{\alpha t}\right)^\frac {1}{\alpha-1} .
$$
Concerning the shock wave solution
$$
u(x,t) = \left\lbrace
\begin{aligned}
&0 &&\text{if}\quad x <st\\
&1 &&\text{if}\quad st <x
\end{aligned}
\right.
$$
with
$$
s = \frac{|1|^\alpha - |0|^\alpha}{1 - 0} = 1
$$
following from the Rankine-Hugoniot condition, one shows that this solution does not satisfy the Lax entropy condition since
$$
\alpha\, 0\, |0|^{\alpha-2} < s < \alpha\, 1\, |1|^{\alpha-2} .
$$
Indeed, the Lax entropy condition requires opposite inequalities.
Finally, the entropy solution is the rarefaction wave.
Best Answer
You are very close, but slightly off, and I think the confusion comes in this variable of integration $r$.
I'll present the exposition that has always made most sense to me. Characteristics $(x(\tau),t(\tau))$ originating from a point $(x(0),t(0)) = (\xi,0)$ are given by \begin{align*} \dot x &= b, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\, x(0) = \xi, \\ \dot t &= 1, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, t(0) = 0, \\ \dot z &= s(x,t), \,\,\, z(0) = u_0(\xi), \end{align*} where $z(\tau;\xi) = u(x(\tau;\xi),t(\tau;\xi))$. Solving these gives \begin{align*} x(\tau;\xi) &= b\tau+\xi, \\ t(\tau;\xi) &= \tau, \\ z(\tau;\xi) &= u_0(\xi) + \int^\tau_{0} s(x(r;\xi),t(r;\xi))dr = u_0(\xi) + \int^\tau_{0} s(br+\xi,r)dr. \end{align*} Now if we can invert the relationships $(x(\tau;\xi),t(\tau;\xi))$ to arrive at $(\tau(x,t),\xi(x,t))$, then we have the solution $u(x,t) = z(\tau(x,t);\xi(x,t))$. Note that the integration variable $r$ will stay! Only $\tau$ and $\xi$ will be replaced. As you've noted $\xi = x-bt$, and we also have $\tau = t$. Thus \begin{align*} u(x,t) &= u_0(x-bt)+\int^t_0 s(br+x-bt,r)dr \\ &= u_0(x-bt)+\int^t_0 s(x-b(t-r),r)dr. \end{align*} We can check that this is indeed a solution. Note that when you differentiate the integral term with respect to $t$, you need to use the Leibniz rule. We see: \begin{align*} u_t(x,t) &= -bu'_0(x-bt) + s(x,t) -b \int^t_0 s_x(x-b(t-r),r)dr,\\ u_x(x,t) &= u'_0(x-bt) + \int^t_0 s_x(x-b(t-r),r)dr. \end{align*} Thus we do indeed have $$u_t(x,t)+bu_x(x,t) = s(x,t),$$ and clearly $u(x,0) = u_0(x)$.