  We perform two-dimensional simulations of Alfven oscillations in magnetars,
modeled as relativistic stars with a dipolar magnetic field. We use the
anelastic approximation to general relativistic magnetohydrodynamics, which
allows for an effective suppression of fluid modes and an accurate description
of Alfven waves. In addition, we compute Alfven oscillation frequencies along
individual magnetic field lines with a semi-analytic approach, employing a
short-wavelength approximation. Our main findings are as follows: a) we confirm
the existence of two families of quasi-periodic oscillations (QPOs), with
harmonics at integer multiples of the fundamental frequency, as was found in
the linear study of Sotani, Kokkotas & Stergioulas; b) the QPOs appearing near
the magnetic axis are split into two groups, depending on their symmetry across
the equatorial plane. The antisymmetric QPOs have only odd integer-multiple
harmonics; c) the continuum obtained with our semi-analytic approach agrees
remarkably well with QPOs obtained via the two-dimensional simulations,
allowing for a clear interpretation of the QPOs as corresponding to turning
points of the continuum. This agreement will allow for a comprehensive study of
Alfven QPOs for a larger number of different models, without the need for
time-consuming simulations. Finally, we construct empirical relations for the
QPO frequencies and compare them to observations of known Soft Gamma Repeaters.
We find that, under the assumptions of our model and if the magnetic field of
magnetars is characterized by a strong dipolar component, and QPOs are produced
near the magnetic pole, then one can place an upper limit to the mean surface
strength of the magnetic field of about 3-8 10^15 G.
