As I said, I haven’t studied Riesz measures for subharmonic functions, so my answer could not be 100% appropriate or canonical.
Concerning your question, the answer is no. The idea is that you are basically multiplying $\phi$ by a characteristic function, so the function you test $\mu_n$ against is not continuous. That is in general a problem for converging measures.
If you take a general succession of Radon measures $\mu_n$ which converges (in measure) to $\mu$, the thesis is not true. For instance, take a radial $\phi$, $\mu$ to be the dirac delta in the origin, $\mu_n$ a standard regularization of it (so, basically, just the standard mollifier), $V=B_1(0)$ and $W=V\cap\{x_1>0\}$. We have
$$\lim_{n\to\infty}\int_W\phi(x)d\mu_n(x)=1/2 \lim_{n\to\infty}\int_V\phi(x)d\mu_n(x)=1/2\int_V\phi(x)d\mu(x),$$
which is in general nonzero, while
$$\int_W\phi(x)d\mu(x)=0.$$
You can realize the measures above with subharmonic functions taking $u=-\Gamma$ and $u_n$ a standard regularization of $u$, where $\Gamma$ is the fundamental solution of the Laplace operator in $\mathbb R^m$, s.t. $\Delta\Gamma=-\delta$ distributionally.
My only doubt is that I don’t know if that $u$ is regular enough for your setting, or it doesn’t count as it is unbounded. If not ok, you could take instead $u=\max(-\Gamma,-c)$ (with $c>0$ large enough for $u$ not to be constant on $V$), that is subharmonic, bounded, and the Laplacian is supported on a spherical surface inside $V$ that is basically the set of points where the gradient of $u$ is discontinuous. Finally, taking $W$ as the complement of that surface should do an analogous trick.
Best Answer
It's better to refer to the divergence theorem (a.k.a. Gauss' theorem), but in two dimensions that follows from Green's theorem: $$ \int_D \nabla \cdot F = \int_D \left( \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} \right) dx \, dy = \{ M= F_x,\ L=-F_y \} \\ = \oint_{\partial D} (-F_y \,dx + F_x \,dy) = \{ n\,ds = (dy, -dx) \} = \oint_{\partial D} F \cdot n \, ds . $$
Then we do a multidimensional integration by parts using the formula $\nabla\cdot(fv) = \nabla f \cdot v + f \, \nabla\cdot v,$ when $f$ is a scalar function and $v$ a vector function: $$ \int_D \phi \, \Delta u = \int_D \nabla\cdot(\phi\,\nabla u) - \int_D \nabla\phi \cdot \nabla u = \underbrace{\oint_{\partial D} (\phi\,\nabla u) \cdot n\,ds}_{=0} - \int_D \nabla\phi \cdot \nabla u . $$ Another integration by parts gives $$ \int_D \nabla\phi \cdot \nabla u = \int_D \nabla\cdot(\nabla\phi \, u) - \int_D \Delta\phi \, u = \underbrace{\oint_{\partial D} (\nabla\phi \, u) \cdot n\,ds}_{=0} - \int_D \Delta\phi \, u . $$
Thus, $$ \int_D \phi \, \Delta u = \int_D \Delta\phi \, u . $$