As discussed in #2452, the integration weights in cylindrical coordinates are implemented identically as those in cartesian coordinates (which is incorrect). The procedure to fix this is outlined in the code itself:

but it's not trivial to implement, particularly with cases at or crossing r=0 (which happens a lot...)

Anecdotally, this doesn't seem to be an issue if R is rather large. However, for smaller cells, the error compounds substantially (and doesn't "easily" disappear with resolution).