I don't think it's so much that this integral vanishes, but rather that it's negative.
This is an application of Green's identity:
$$\begin{align}
\int_\Omega v\cdot\Delta v\,dx & = -\int_\Omega |\nabla v|^2dx + \int_{\partial\Omega}v\,\nabla v\cdot n\,dx \\
& = -\int_\Omega|\nabla v|^2dx \\
\end{align}$$
The boundary integral vanishes because we assumed that $v = u_1 - u_2$ is the difference of two hypothetical solutions satisfying the same Dirichlet boundary conditions.
So to fill out the chain of inequalities from the original post just a bit:
$$\begin{align}
\frac{1}{2}\frac{d}{dt}\int_\Omega v^2dx & = \int_\Omega v\cdot\partial_tv\, dx \\
& = \int_\Omega v\cdot(\partial_t v - \Delta v + \Delta v)dx \\
& \le \int_\Omega v\cdot\Delta v\,dx + \int_\Omega|v||\partial_t v - \Delta v|dx \\
& = -\int_\Omega|\nabla v|^2dx + \int_\Omega|v||\partial_t v - \Delta v|dx \\
& \le \int_\Omega|v||\partial_tv - \Delta v|dx
\end{align}$$
and from there you apply the rest of the specific knowledge from the problem, like the fact that the nonlinear term is cosine which is a Lipschitz function.
To use a little more fancy terminology, the Laplace operator acting on function spaces where any of the usual boundary conditions apply (Dirichlet, Neumann, or Robin) is symmetric and negative-definite.
This is also true of other divergence-form elliptic operators; it's tremendously useful for proving stability properties about the solution and for numerical analysis.
