How do I solve a 3D poisson equation with mixed neumann and periodic boundary conditions numerically? The PDE is being solved over a cube. Four of the faces are periodic, and one set of opposing faces have no-flow Neumann boundary conditions:
$\nabla^2 u = f(x,y,z)$
$\frac{\partial u}{\partial z} = 0$ at $z = 0$ and $z = h$
I am able to solve a similar problem with full Neumann boundary conditions (assuming the compatibility condition for $f(x,y,z)$ is satisfied) by arbitrarily selecting one node and setting it to $0$ as suggested here How to numerically solve the Poisson equation given Neumann boundary conditions? 
So, based on this I have 2 questions:


*

*Would this method also work for mixed periodic/Neumann boundary conditions and if so why? I can't seem to find an answer. I believe periodic boundary conditions also have a family of solutions valid up to an additive constant $C$.

*Surely the choice of which node you choose to set to $0$ has an affect on the final solution. Therefore, are there any rules of thumb on how to go about this? Or am I understanding the problem wrong?
I think this is a generally applicable problem, so any answers would be appreciated.
 A: $\nabla^2$ is separable in space:
$$\nabla^2 = \nabla^2_x + \nabla^2_y + \nabla^2_z$$
where $ \nabla^2_x + \nabla^2_y$ are with the periodic condition, and $\nabla_z^2$ be with the replicate condition. Notice that for discrete 1D Laplacian operator $[1, -2, 1]$, replicate and symmetric boundary conditions are the same. So along $z$ you can safely assume symmetric boundary condition.
Given these information, you can do Fourier transform along $x$ and $y$ to get the inverse of $\nabla^2_x + \nabla^2_y$, and discrete cosine transforms along $z$ to get the inverse of $\nabla^2_z$.
Denote $\mathcal{F}_x$ and $\mathcal{F}_y$ as Fourier transforms along $x$ and $y$, $\mathcal{D}_z$ as discrete cosine transform along $z$, your solution can be obtained as:
$$u = \mathcal{D}_z^{-1}\left\{\mathcal{F}^{-1}_y\left\{\mathcal{F}^{-1}_x \left\{ \frac{\mathcal{D}_z\{\mathcal{F}_y\{\mathcal{F}_x\{f\}\}\}}{\mathcal{D}_z\{\mathcal{F}_y\{\mathcal{F}_x\{\nabla^2\}\}\}} \right\}\right\}\right\}$$
As for your questions:


*

*I think it is Yes.

*For this (direct) approach, you have direct point-to-point division in the spectral domain. To avoid singularity at the DC term (the mean value), you can set any non-zero number. It only affects the mean of the solution.

