  In axisymmetric fusion reactors, the equilibrium magnetic configuration can
be expressed in terms of the solution to a semi-linear elliptic equation known
as the Grad-Shafranov equation, the solution of which determines the poloidal
component of the magnetic field. When the geometry of the confinement region is
known, the problem becomes an interior Dirichlet boundary value problem. We
propose a high order solver based on the Hybridizable Discontinuous Galerkin
method. The resulting algorithm (1) provides high order of convergence for the
flux function and its gradient, (2) incorporates a novel method for handling
piecewise smooth geometries by extension from polygonal meshes, (3) can handle
geometries with non-smooth boundaries and x-points, and (4) deals with the
semi-linearity through an accelerated two-grid fixed-point iteration. The
effectiveness of the algorithm is verified with computations for cases where
analytic solutions are known on configurations similar to those of actual
devices (ITER with single null and double null divertor, NSTX, ASDEX upgrade,
and Field Reversed Configurations).
