  We present an adaptive multilevel Monte Carlo (AMLMC) algorithm for
approximating deterministic, real-valued, bounded linear functionals that
depend on the solution of a linear elliptic PDE with a lognormal diffusivity
coefficient and geometric singularities in bounded domains of $\mathbb{R}^d$.
  Our AMLMC algorithm is built on the results of the weak convergence rates in
the work [Moon et al., BIT Numer. Math., 46 (2006), 367-407] for an adaptive
algorithm using isoparametric d-linear quadrilateral finite element
approximations and the dual weighted residual error representation in a
deterministic setting. Designed to suit the geometric nature of the
singularities in the solution, our AMLMC algorithm uses a sequence of
deterministic, non-uniform auxiliary meshes as a building block. The
deterministic adaptive algorithm generates these meshes, corresponding to a
geometrically decreasing sequence of tolerances. For a given realization of the
diffusivity coefficient and accuracy level, AMLMC constructs its approximate
sample using the first mesh in the hierarchy that satisfies the corresponding
bias accuracy constraint.
  This adaptive approach is particularly useful for the lognormal case treated
here, which lacks uniform coercivity and thus produces functional outputs that
vary over orders of magnitude when sampled.
  We discuss iterative solvers and compare their efficiency with direct ones.
To reduce computational work, we propose a stopping criterion for the iterative
solver with respect to the quantity of interest, the realization of the
diffusivity coefficient, and the desired level of AMLMC approximation.
  From the numerical experiments, based on a Fourier expansion of the
coefficient field, we observe improvements in efficiency compared with both
standard Monte Carlo and standard MLMC for a problem with a singularity similar
to that at the tip of a slit modeling a crack.
