Weak solutions to $\Delta u=f$ are in $W^{2,2}$ I  believe  the following statement is true.
Let $\Omega$ be a smoothly, bounded domain in $\mathbb{R}^{n}$.
The statement:
Let $u\in H^{1}(\Omega)$ so that
there exists $f\in L^{2}(\Omega)  \;s.t.\int_{\Omega}\nabla u\nabla \varphi=\int_{\Omega} f\varphi, \forall \varphi\in H^{1}(\Omega)$.
Then $u\in H^{2}(\Omega)$.

I have searched in the book by Evans and Brezis but not so certain. Could anyone provide a reference for that? It can be seen in Evans's book that $u\in H^{2}_{loc}(\Omega)$.

Thanks so much.

Please check a general elliptic regularity result ($L^p$ version) on: 

Dauge, Monique. Elliptic boundary value problems on corner domains: smoothness and asymptotics of solutions, 1988

Theorem 20.10 (together with the explanation of notation above the theorem), essentially we have the following regularity result:
$$
\|u\|_{H^2} \leq C \|\Delta u\|_{L^2}.
$$
i.e., in your case, $\Delta u = f$ which is implied the fundamental theorem of calculus of variation. For a more specific version in $L^2$ sense, please refer to Grisvard's book Elliptic Problems in Nonsmooth Domains, $\S 2.3.3$, where the whole chapter 2 deals with weak solutions.
BTW: please be careful with the test function you choose, if your test function is in $H^1(\Omega)$, not in $H^1_0(\Omega)$, this is a Neumann problem and your $u$ and $f$ should satisfy certain compatibility condition. 
