First of all, as the function 1/tan(x) is singular at at multiple points in your domain, it does not suprise me that you get a nan answer to the PDE. (This means that when cot is interpolated onto a finite element space, one of it values are nan).
This can at least be partially circumvented by letting your mesh not include the zero node, like:
lower=1e-5
upper=2*np.pi
mesh=IntervalMesh(1000,lower,upper)
and evaluate the solution at u0(lower).