First, what does it mean to say
NeumannValue[r, z == 0]
taken into account that r is an independent variable?
Then there are not enough conditions to avoid randomness in a way. To see that take manual's example
NDSolveValue[{\!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[x, y]\)\) ==
NeumannValue[1., x >= 0.35],
DirichletCondition[u[x, y] == 0., x <= -0.3]}, u, {x, y} \[Element]
Disk[]]
Plot3D[%[x, y], {x, y} \[Element] Disk[]]
and kick the Dirichlet condition out. The result is (in 2D)
In[9]:= NDSolveValue[{\!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[x, y]\)\) ==
NeumannValue[1., x >= -1.]}, u, {x, y} \[Element] Disk[]]
During evaluation of In[9]:= NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified; the result may be off by a constant value. >>
Out[9]= InterpolatingFunction[{{-1., 1.}, {-1., 1.}}, <>]
and this NDSolveValue::femibcnd is pretty reasonable. In 3D the question is: Where are the conditions for the shank?