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.
The regularity theorem which states that if $\color{blue}{f\in L^2(\Omega)}$ and $$ \color{blue}{\int_{\Omega}\nabla u\nabla \varphi\,\mathrm{d}x=\int_{\Omega} f\varphi\,\mathrm{d}x\quad \forall \varphi\in H^{1}(\Omega)}\implies\color{blue}{u\in H^2(\Omega)} \tag{PR}\label{pr} $$ is true if and only if some regularity conditions on the boundary of the domain $\partial \Omega$ and on the boundary value of the solution $u$ are assumed (be it a Dirichlet, Neumann or Robin boundary value problem). Otherwise, as @ArcticChar states, you can only expect $u\in H_\mathrm{loc}^2(\Omega)$.
Precisely
The first step of Mikhailov is to prove the theorem with homogeneous boundary conditions:
Theorem 4 ([2], p. 217). If $f\in H^k(\Omega)$ and $\partial\Omega\in C^{k+2}$ for certain $k\ge 0$, then the generalized solutions $u(x)$ of the first and second boundary-value problems (respectively called Dirichlet and Neumann problems) with homogeneous boundary condition for the Poisson equation belong to $H^{k+2}(\Omega)$ and satisfy (in the case of the second boundary problem it is assumed that $\int_\Omega u \mathrm{d}x=0$) the inequality $$ \Vert u\Vert_{ H^{k+2}(\Omega)}\le C\Vert f\Vert_{ H^{k}(\Omega)} $$ where the constant $C>0$ does not depend on f.
Then he uses the regularity theorem 4 above to the case of non homogeneous boundary conditions ([2], p. 226) by reducing to homogeneous ones. He explicitly does this only for the Dirichlet problem as we do, but the same method nevertheless works for the Neumann problem, as shown in this answer. Consider a solution $u(x)$ of the problem Dirichlet problem above and a function $\Phi\in H^{k+2}(\Omega)$ whose trace on $\partial\Omega$ is $g\in H^{k+1/2}(\partial\Omega)$, i.e. $$ \Phi=g\;\text{ on }\;\partial\Omega\label{1}\tag{1} $$ Define $w=u-\Phi$: then, for every $\varphi\in H^{1}_0(\Omega)$, $$\label{GeneralizedDirichlet}\tag{GN} \begin{split} \int_{\Omega} \nabla w \cdot \nabla \varphi \, \mathrm{d}x&=\int_{\Omega} \nabla u \cdot \nabla \varphi \, \mathrm{d}x - \int_{\Omega} \nabla \Phi \cdot \nabla \varphi \, \mathrm{d}x \\ &= \int_\Omega f \varphi \, \mathrm{d}x + - \int_{\Omega} \nabla \Phi \cdot \nabla \varphi \, \mathrm{d}x \\ &=\int_\Omega f \varphi \, \mathrm{d}x - \int_{\Omega} \nabla\cdot(\varphi\nabla\Phi) \, \mathrm{d}x + \int_{\Omega} \varphi\Delta\Phi \, \mathrm{d}x \\ &=\int_\Omega \big[\,f + \Delta\Phi\big]\varphi\, \mathrm{d}x , \end{split}\label{2}\tag{2} $$ thus we get $$ \int_{\Omega} \nabla w \cdot \nabla \varphi \, \mathrm{d}x=\int_\Omega f_1 \varphi \, \mathrm{d}x $$ where $f_1=f+\Delta\Phi\in H^k(\Omega)$. This means that $w$ solves the homogeneous Dirichlet problem with homogeneous boundary conditions (i.e. $g\equiv 0$), and by applying theorem 4 we have $$ w=u-\Phi\in H^{k+2}(\Omega) \iff u=w+\Phi\in H^{k+2}(\Omega) $$
Notes
[1] P. Grisvard (1985), Elliptic problems in nonsmooth domains, Monographs and Studies in Mathematics, 24. Pitman Advanced Publishing Program. Boston-London-Melbourne: Pitman Publishing Inc., pp. XIV+410, MR0775683, Zbl 0695.35060.
[2] V. P. Mikhailov (1978), Partial differential equations, Translated from the Russian by P.C. Sinha. Revised from the 1976 Russian ed., Moscow: Mir Publishers, p. 396 MR0601389, Zbl 0388.3500.