Solution (source code)

= Solution

For $\phi\in C_0^\infty(D)$, set
$$
f_\phi=-2\pi\Delta^{-1}\phi
=\int_DG_D(\mathord\cdot,y)\phi(y)\,dy
\in H_0^1(D)
$$
and define the distributional pairing by
$$
(h,\phi):=(h,f_\phi)_\nabla.
$$
It is centered Gaussian by the definition of the GFF. Integration by parts gives
$$
\begin{aligned}
\operatorname{Var}(h,\phi)
&=\|f_\phi\|_\nabla^2\\
&=\frac1{2\pi}\int_Df_\phi(-\Delta f_\phi)\,dx\\
&=\int_Df_\phi(x)\phi(x)\,dx\\
&=\iint_{D\times D}\phi(x)G_D(x,y)\phi(y)\,dx\,dy.
\end{aligned}
$$
This is the <Test-function pairing with a Gaussian free field>.