How to compute the gradient of pde constrained optimal control problem with Neuman boundary control

This problem is not unique if you only have Neumann conditions.
See: Issues Solving Pure Neumann Problems in dolfinx - #3 by jkrokowski

I also dont think your Adjoint problem is correct, as you have changed your homogeneous Neumann conditions into Dirichlet conditions.