$\nabla f = h$ is intended in the sense of distributions, i.e. $\int f \mathrm{div}\, \phi = -\int h \phi$ for all smooth compactly supported $\phi$.

@Drazick problem solved in the meantime?

Can you specify your boundary condition, and your assumptions on the initial data?

@delimit No, $C_c^\infty$ is not dense in $L^\infty(L^\infty)$. But then you can do the approximation with weaker information than that. For instance, it is enough that the approximating sequence $u_m$ (a) converges pointwise, (b) has a uniform bound in $L^\infty(L^\infty)$, (c) has $\nabla u_m$ converging in $L^2(L^2)$, and (d) has $\partial_t u_m$ converging in $L^2(H^{-1})$. Typically approximation by convolution will do that for you.

@YemonChoi, thinking about the injection as a quotient map makes things much clearer. Thanks!

@SeanEberhard, that is very useful! As I understand it, $C_0$ is not dense in $C_b$, and therefore knowing $\xi\in C_b'$ on $C_0$ is not enough to determine $\xi$.

Do you know of a good reference for smoothing properties of vectorial Lp spaces?

Shouldn't the second bracket read $\langle v,f'(u_\epsilon)u'_\epsilon \rangle$?

If the brackets are in space, then it seems the full time derivative is missing ..

Do you mean the brackets to be in space or space-time? I'm trying to understand the meaning of the displayed equation ...

I just realized that your case is slightly better than the standard $L^1$-right-hand-side case, since you have bounded $\|\nabla \Psi_n\|_2$ by construction. By the Sobolev inequality that implies that $\Psi_n$ is bounded in $L^{d/(d-2)}$ and you can extract a weakly converging subsequence. Then you can pass to the limit in your equality $\Psi_n = G*u_n$ after integrating against a test function.