Show that there exists a unique weak solution $(u, p)$ such that $-\Delta u + p = f$ Let $\Omega\subseteq\mathbb{R}^n$ be a bounded domain and $f\in L^2(\Omega)$.
I would like to show that there exists a unique pair $(u, p)$ with $u \in H_0^1(\Omega)$ and $p\in \mathbb{R}$ such that
$$
\int_\Omega \nabla u\cdot \nabla v + p\int_\Omega v = \int_\Omega fv \qquad\text{and}\quad\quad \int_\Omega u = 0
$$
for all $v\in C_c^\infty(\Omega)$.
If we fix $p$ and do not require that $\int_\Omega u = 0$, then this problem is a straightforward problem that can be solved using well-known methods in PDE.
My idea was therefore to show that we must have $p = \int_\Omega f$. This would be easy to show by setting $v=1$ if it were not for the requirement $v\in C_c^\infty(\Omega)$. Then, for this choice of $p$ I can find a solution $u\in H_0^1(\Omega)$ such that
$$
\int_\Omega \nabla u\cdot \nabla v + p\int_\Omega v = \int_\Omega fv
$$
for all $v\in C_c^\infty(\Omega)$.
I was then hoping to conclude that the condition $\int_\Omega u=0$ must also be satisfied.
One possible approach is to consider $H_0^1(\Omega)\times \mathbb{R}$ as a Hilbert space and apply the Lax-Milgram theorem to an appropriate choice of bilinear map and bounded linear functional. More specifically, I consider the bilinear form
$$
B((u,p), (v,q)) = \int_\Omega fv + p\int_\Omega v - q\int_\Omega u + pq
$$
and the linear function
$$
L((v,q)) = \int_\Omega f(v-q).
$$
Here, $u,v\in H_0^1(\Omega)$ and $p,q\in \mathbb{R}$. I can then conclude that there exists a unique pair $(u, p)$ satisfying
$$
\int_\Omega fv + p\int_\Omega v - q\int_\Omega u + pq
=\int_\Omega f(v-q)
$$
for all $v\in H_0^1(\Omega)$ and $q\in \mathbb{R}$.  In particular, setting $q=0$ we see that
$$
\int_\Omega fv + p\int_\Omega v
=\int_\Omega fv
$$
as desired. Furthermore, ranging $q\in \mathbb{R}$ I can conclude that
$$
\int_\Omega u = p - \int_\Omega f.
$$
However, I cannot show that $\int_\Omega u = 0$.
 A: Probably the easiest way to prove it is to see that it is equivalent to optimizing the Dirichlet functional with the constraint $\int_\Omega u\, dx = 0$, prove that it has a minimum using a minimizing sequence and then applying Lagrange theorem to obtain a solution of your problem ($p$ would be the Lagrange multiplier associated to the restriction), but in this case the constraint is easy enough that we can do it more directly.
Let us consider the subspace
$$ V = \left \{u \in H_0^1(\Omega)| \int_\Omega u = 0\right\},$$
which is a closed subspace of $H_0^1(\Omega)$, hence is it Banach, and if we define the usual bi-linear Laplacian form
$$ B(u,v) = \int_\Omega \nabla u \nabla v\, dx,$$
and the linear form
$$ \ell(v) = \int_\Omega f v\, dx,$$
They still satisfy all the hypotesis of Lax-Milgram theorem (or just Riesz representation theorem) because they satisfy them in the bigger space $H_0^1(\Omega)$. Therefore we know that there is a unique $u\in V$ such that
$$ \int_\Omega \nabla u \nabla v \, dx = \int_\Omega fv \, dx, \quad \forall v\in V.$$
Now we need to find out how this weak formulation looks when $v \in H_0^1(\Omega)$, but as you already noticed, we cannot simply substract the average from $v$ because $1\not\in H_0^1(\Omega)$. Instead we consider a fixed function $\phi \in C_0^\infty(\Omega)$ such that $\int_\Omega \phi \, dx = 1$ depending only on the domain $\Omega$.
Then, for every $v\in H_0^1(\Omega)$, we have that $v - 
\phi \bar v \in V$ where
$$\bar v = \frac{1}{|\Omega|} \int_\Omega v \, dx,$$
is the average of $v$.
Hence we obtain
$$ \int_\Omega \nabla u \nabla v \, dx  - \bar v \int_\Omega \nabla u \nabla \phi \, dx = \int_\Omega fv \, dx - \bar v\int_\Omega f\phi \, dx , \quad \forall v\in H_0^1(\Omega).$$
which can be writen as
$$ \int_\Omega \nabla u \nabla v \, dx  + p\int_\Omega v \, dx = \int_\Omega fv \, dx, \quad \forall v\in H_0^1(\Omega),$$
with
$$p = \frac{1}{|\Omega|}\int_\Omega f\phi - \nabla u\nabla \phi \, dx$$
which is equivalent to your weak formulation due to the density of $C_0^\infty(\Omega)$ in $H_0^1(\Omega)$.
For the uniqueness, Lax-Milgram (or Riesz) gives us uniquenes of $u$ and we saw that $u$ completely termins $p$ so it is also unique (here is important that even though $\phi$ is somewhat arbitrary, it only depends on the domain).
Remarks

*

*If the fact that the expression for $p$ has the function $\phi$ on it (even though it doesn't depend on it), at least if the domain is regular enough we can use the regularity the theorem: We have that $f - p \in L^2(\Omega)$ so $u\in H^2(\Omega)$, and then we have that $u$ satisfies
$$ -\Delta u + p = f, $$
not only as distributions but on $L^2(\Omega)$, and after integrating both sides, and using the divergence theorem,
$$p = \bar f - \frac{1}{|\Omega|} \int_{\partial \Omega} \frac{\partial u}{\partial n} \, dS. $$


*If the boundary is not regular enough, I don't know if you can find a expression for $p$ that doesn't involve $\phi$, but you can prove that doesn't depend on it. If $\psi$ is another of such functions, then $\phi - \psi \in V$
$$ \int_\Omega f(\phi - \psi) - \nabla u\nabla (\phi - \psi) \, dx = 0,$$
from where it follows
$$ \int_\Omega f\phi- \nabla u\nabla \phi \, dx = \int_\Omega f\psi- \nabla u\nabla \psi \, dx. $$
A: Uniqueness only
Suppose $(u, p)$ and $(u', p')$ are two solutions and set $w = u - u'$, $q = p - p'$. Then $\int_\Omega\nabla w \cdot \nabla v + q \int_\Omega v = 0$ for all test functions $v$. With $v = w$ you have $\int_\Omega \vert\nabla w\vert^2 + q \int_\Omega w = 0$. But the latter integral is zero since $\int_\Omega u = \int_\Omega u' = 0$, so $\nabla w$ = 0, that is $w = \text{const}$. With $\int_\Omega w = 0$ you conclude that the constant is zero, thus $u = u'$. Call this $u_0$.
Now suppose $(u_0, p)$ and $(u_0, p')$ are two solutions. It follows
\begin{align*}
\int_\Omega\, \nabla u_0 \cdot \nabla v + p\int_\Omega\ v &= \int_\Omega\, fv\\
\int_\Omega\, \nabla u_0 \cdot \nabla v + p'\int_\Omega\ v &= \int_\Omega\, fv
\end{align*}
Subtracting the second from the first gives $(p - p')\int_\Omega v = 0$ for all $v$.
Sleight of hand: You actually need to use a density/compactness argument with approximations of $u$, $u'$, $u_0$ arbitrarily close to elements of $C^\infty_c(\Omega)$.
