  The equilibrium configuration of a plasma in an axially symmetric reactor is
described mathematically by a free boundary problem associated with the
celebrated Grad--Shafranov equation. The presence of uncertainty in the model
parameters introduces the need to quantify the variability in the predictions.
This is often done by computing a large number of model solutions on a
computational grid for an ensemble of parameter values and then obtaining
estimates for the statistical properties of solutions. In this study, we
explore the savings that can be obtained using multilevel Monte Carlo methods,
which reduce costs by performing the bulk of the computations on a sequence of
spatial grids that are coarser than the one that would typically be used for a
simple Monte Carlo simulation. We examine this approach using both a set of
uniformly refined grids and a set of adaptively refined grids guided by a
discrete error estimator. Numerical experiments show that multilevel methods
dramatically reduce the cost of simulation, with cost reductions typically on
the order of 60 or more and possibly as large as 200. Adaptive gridding results
in more accurate computation of geometric quantities such as x-points
associated with the model.
