The gyrokinetic Poisson equation arises as a subproblem of Tokamak fusion reactor simulations. It is often posed on disk-like cross sections of the Tokamak that are represented in generalized polar coordinates. On the resulting curvilinear anisotropic meshes, we discretize the differential equation by finite differences or low order finite elements. Using an implicit extrapolation technique similar to multigrid $$\tau $$
τ
-extrapolation, the approximation order can be increased. This technique can be naturally integrated in a matrix-free geometric multigrid algorithm. Special smoothers are developed to deal with the mesh anisotropy arising from the curvilinear coordinate system and mesh grading.