Development of a new zonal flow equation solver by diagonalisation and its application in non-circular cross-section tokamak plasmas. A toroidal gyrokinetic full-f code GKNET (GyroKinetic Numerical Experimental Tokamak) with field solver in real space has been developed recently to simulate micro-turbulence dynamics in the circular cross-section tokamak plasmas (Obrejan et al., 2015). In this work, we introduce a new high accuracy Zonal Flow (ZF) equation solver which makes use of a parametrisation of the D-shaped magnetic flux surfaces to diagonalise the ZF equation. In addition to being more rigorous near the magnetic axis of the poloidal plane compared to methods based on local approximations, the ZF solver here allows to properly take into account the shape of magnetic flux surfaces independently of the coordinate system used in the rest of the code. The upgraded GKNET code is applied to study the collisionless damping of the Geodesic Acoustic Modes (GAMs) in elliptic and both positive and negative D-shaped magnetic configurations. We found that in addition to the influence of elongation, triangularity is effective in increasing the damping rate of GAMs, independently of the sign of the triangularity.