  We are given a uniformly elliptic coefficient field that we regard as a
realization of a stationary and finite-range (say, range unity) ensemble of
coefficient fields. Given a (deterministic) right-hand-side supported in a ball
of size $\ell\gg 1$ and of vanishing average, we are interested in an algorithm
to compute the (gradient of the) solution near the origin, just using the
knowledge of the (given realization of the) coefficient field in some large box
of size $L\gg\ell$. More precisely, we are interested in the most seamless
(artificial) boundary condition on the boundary of the computational domain of
size $L$.
  Motivated by the recently introduced multipole expansion in random media, we
propose an algorithm. We rigorously establish an error estimate (on the level
of the gradient) in terms of $L\gg\ell\gg 1$, using recent results in
quantitative stochastic homogenization. More precisely, our error estimate has
an a priori and an a posteriori aspect: With a priori overwhelming probability,
the (random) prefactor can be bounded by a constant that is computable without
much further effort, on the basis of the given realization in the box of size
$L$.
  We also rigorously establish that the order of the error estimate in both $L$
and $\ell$ is optimal, where in this paper we focus on the case of $d=2$. This
amounts to a lower bound on the variance of the quantity of interest when
conditioned on the coefficients inside the computational domain, and relies on
the deterministic insight that a sensitivity analysis wrt a defect commutes
with (stochastic) homogenization. Finally, we carry out numerical experiments
that show that this optimal convergence rate already sets in at only moderately
large $L$, and that more naive boundary conditions perform worse both in terms
of rate and prefactor.
