I am playing around with implementing an elliptic divergence cleaning scheme for induction equation in MHD. My grid is staggered, i.e. the magnetic field components $B_x, B_y, B_z$ are at cell faces while the scalar $\psi$ (of the poisson equation $(\nabla^2 \psi = -\nabla \cdot B)$ for mag field correction) is computed at cell centers. This way, the gradient of $\psi$ is computed at cell faces, which when added to the old $B^o_x$ value, yields a new, updated $B^n_x$, as,
$B^n_x = B^o_x + \nabla_x \psi$
where $B^n_x$ is updated to be divergence free.
Now my problem is that on the boundary node, I cannot correct for $B_x$, since I cannot compute $\nabla_x \psi$ at the boundary. Any inputs to fix this? The boundary condition on $B_x$ is Neumann, $\partial_x B_x |_{wall} = 0$.