differential equations - Accuracy of Grad-Shafranov PDE's NDSolveValue over implicit region

Context

I need to solve for the toroidal flux of the magnetic field above an accretion disc.

For this purpose, I define the region over which the flux is non zero,

a= 2;
f[R_] = a (1 - R^2)^(1/2 + 1/100);
Ω = ImplicitRegion[z <= f[R], {{R, 0, 1}, {z, 0, f[0]}}];

Then I solve for the force-free Grad-Shafranov equation:

eqn0 =R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] - D[P[R, z], R] == -R/2;

as follows

P0 = NDSolveValue[{eqn0, DirichletCondition[P[R, z] == 0, R == 0], 
    DirichletCondition[P[R, z] == 0, z == f[R]]},  P, {R, z} ∈ Ω];

I then plot the resulting (normalized) flux map:

np = NMaximize[P0[R, z], {R, z} ∈ Ω][[1]];
ContourPlot[ P0[R, z]/np, {R, z} ∈ Ω,PlotPoints -> 50, 
ImageSize -> Small, AspectRatio -> Sqrt[f[0]]]

and look at the value of the flux on the boundary

Plot[P0[R, z]/np /. z -> f[R], {R, 0, 1}, PlotRange -> All]

it seems to satisfy the boundary condition.

If I now look at the pressure above the cap:

grad2 = Grad[P0[R, z]/np, {R, z}] // (#.#/R^2 &);
Plot[grad2 /. z -> f[R], {R, 0, 1}]

It is not smooth enough…

On top of that

If I decide to extend the height of the column over which the flux is defined, to say

 a=15; 

Then the accuracy of he map deteriorates considerably for the map

and even more for the pressure:

QUESTION(S)

How can improve the accuracy of the solution found by NDSolveValue?

I understand that there are options such as PrecisionGoal or Method, but I guess I am trying to ask a more general question:

More generally, what is the best learning strategy within mathematica to be able to find such improvement? (a.k.a how not to get lost in the documentation?). The idea being, next time I can figure this myself :-)

Mathematica gives the following warning, which is undoubtedly a hint

NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

but I do not know how to follow it up.

PS: If, instead of

 f[R_] = a (1 - R^2)^(1/2 + 1/100);

I have

f[R_] = a (1 - R^2)^(1/2);

the integrator also fails miserably, which is rather odd.