  We present a novel solver technique for the anisotropic heat flux equation,
aimed at the high level of anisotropy seen in magnetic confinement fusion
plasmas. Such problems pose two major challenges: (i) discretization accuracy
and (ii) efficient implicit linear solvers. We simultaneously address each of
these challenges by constructing a new finite element discretization with
excellent accuracy properties, tailored to a novel solver approach based on
algebraic multigrid (AMG) methods designed for advective operators. We pose the
problem in a mixed formulation, introducing the directional temperature
gradient as an auxiliary variable. The temperature and auxiliary fields are
discretized in a scalar discontinuous Galerkin space with upwinding principles
used for discretizations of advection. We demonstrate the proposed
discretization's superior accuracy over other discretizations of anisotropic
heat flux, achieving error $1000\times$ smaller for anisotropy ratio of $10^9$,
for $closed$ $field$ $lines$. The block matrix system is reordered and solved
in an approach where the two advection operators are inverted using AMG solvers
based on approximate ideal restriction (AIR), which is particularly efficient
for upwind discontinuous Galerkin discretizations of advection. To ensure that
the advection operators are non-singular, in this paper we restrict ourselves
to considering open (acyclic) magnetic field lines for the linear solvers. We
demonstrate fast convergence of the proposed iterative solver in highly
anisotropic regimes where other diffusion-based AMG methods fail.
