It looks like specifying the boundary conditions with NeumannValue and DirichletCondition fixes the problem
NDSolve[{Laplacian[u[r], {r, th, z}, "Cylindrical"] + dp == 
   NeumannValue[0, r == 0], 
  DirichletCondition[u[r] == 0, r == 1]}, u, {r, 0, 1}]
But Laplacian has to be used because this didn't work:
NDSolveValue[{1/r D[r u'[r]] + dp == NeumannValue[0, r == 0], 
   DirichletCondition[u[r] == 0, r == 1]}, u, {r, 0, 1}]
So, is there another way around this? My other equations are something like 
1/r D[r U'[r]]/r - dP == 0
1/r D[r k'[r],r] + U'[r]^2 - e[r] == 0
1/r D[r e'[r],r] + e[r]/k[r] U'[r]^2 - c2 re e[r]^2/k[r] ==  0
Where k[r]=0 at r=0 but e[r] goes to zero faster