Suppose $u \in L^2(0,T;L^2)$, $u_t \in L^2(0,T;H^{-1})$ and $f \in L^\infty((0,T)\times\Omega)$. I have the weak form $$\langle u_t, \varphi \rangle_{H^{-1}, H^1} + \int_\Omega\nabla u \nabla \varphi = \int_\Omega f\varphi$$ for all $\varphi \in L^2(0,T;H^1)$. Heuristically, let us pick $\varphi = u_t$: $$|u_t|_{L^2}^2 + \frac{1}{2}\frac{d}{dt}\int_\Omega|\nabla u|^2\leq |f|_{L^2}|u_t|_{L^2}$$ Now using Young's inequality on the RHS, we can cancel out the $|u_t|_{L^2}$ terms and we get $$\frac{d}{dt}\int_\Omega|\nabla u|^2\leq C_f.$$ This gives us a bound on $|\nabla u|_{L^\infty(0,T;L^2)}$.
Obviously this calculation is invalid because we can't pick $\varphi = u_t$ since $u_t$ is not $H^1$ in space, and also we used $\frac{1}{2}\frac{d}{dt}|\nabla u|^2 = \nabla u \cdot \nabla u_t$ which is invalid for the same reason.
How do I do make this calculation rigorous? I was thinking of using approximation by density, but I don't think that works since what equation do the approximating sequence elements satisfy?