# Multilevel Monte Carlo methods for the Grad-Shafranov free boundary problem

Howard C. Elman<EMAIL_ADDRESS>Department of Computer Science and Institute
for Advanced Computer Studies, University of Maryland, College Park.  Jiaxing
Liang<EMAIL_ADDRESS>Applied Mathematics & Statistics and Scientific
Computing Program, University of Maryland, College Park.  Tonatiuh Sánchez-
Vizuet<EMAIL_ADDRESS>Department of Mathematics, The University of
Arizona.

###### Abstract

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.

###### keywords:

Multilevel Monte Carlo Finite-Element , Uncertainty Quantification , Free
Boundary Grad-Shafranov problem , Adaptive Finite Element Discretization.

###### MSC:

[2020] 65Z05 , 65C05 , 62P35 , 35R35 , 35R60

theorem

##  1 Introduction

Monte Carlo (MC) techniques are one of the most common strategies for dealing
with the quantitative assessment of the accuracy of numerical simulations of
physical models with uncertainties. The idea behind these methods is to obtain
a large number of samples (typically by numerically solving the associated
model) for random realizations of the uncertain parameters, and use these data
to gather statistical information about the quantity of interest. However,
when the model involves the solution of partial differential equations, the
computational effort related to the collection of the data points required can
easily become unmanageable. To overcome this difficulty methods like
polynomial chaos expansions [46], stochastic Galerkin [18], and stochastic
collocation [3] have been used to handle uncertainties associated with a small
number of parameters. These techniques, however, often require the development
of specialized numerical solvers or rely on the smooth dependence of the model
with respect to the parameter values. The multilevel Monte Carlo (MLMC) method
was developed [8, 19, 41, 44] as an efficient alternative that does not
require additional smoothness assumptions and can take advantage of an
existing numerical solver. Given a target numerical grid (i.e. a grid whose
resolution is considered sufficiently fine) MLMC improves the efficiency of
the sampling step by offsetting the bulk of the numerical computations to a
sequence of coarser grids where the numerical solution is cheaper.

In the particular context of free-boundary Grad-Shafranov computations subject
to parameter uncertainty, the authors have shown that the computational cost
can be reduced manyfold by employing a strategy based on stochastic
collocation [12]; however, due to the latent possibility of plasma-wall
contacts, the smoothness of the mapping between coil currents and equilibria
cannot be guaranteed. In this paper, our goal is to overcome this difficulty
by approximating the expectation of the equilibrium configuration using a
multilevel Monte Carlo Finite-Element (MLMC-FE) approach. We will consider two
MLMC-FE algorithms: a classical strategy based on uniformly refined meshes and
a variation based on meshes refined adaptively using an a posteriori error
estimator. As we shall see, both of these approaches greatly reduce
computational costs, with the adaptive strategy being somewhat more effective
in computing the approximation of the expectation of geometric properties of
the equilibrium configuration.

An outline of the paper is as follows. In Section 2, we briefly recall the
Grad-Shafranov free boundary problem. Section 3 is devoted to the introduction
of the Monte Carlo and multilevel Monte Carlo Finite-Element methods. The
section concludes with the introduction of an algorithm to compute the optimal
number of samples required at each discretization level. In Section 4 we
present numerical experiments comparing the effectiveness of these strategies
with a Monte Carlo strategy based on a single mesh. Concluding remarks are
presented in the final Section 5. For completeness, the technical
mathematical, and algorithmic aspects of the problem and the methods discussed
are included as an appendix.

##  2 The Grad-Shafranov free boundary problem

###  2.1 The deterministic problem

Figure 1: The plasma confinement region $\Omega_{p}$ is enclosed by the violet
line. The rectangles represent the external coils $C_{k}$; the grey curved
line represents the exterior wall of the vacuum chamber.

In a cylindrically symmetric magnetic confinement device, with coordinates
$(r,z,\varphi)$, the mathematical expression of the equilibrium condition
between the hydrostatic and magnetic forces acting on the plasma results in
the celebrated Grad-Shafranov equation [21, 36, 42]. This nonlinear elliptic
equation relates the poloidal flux function $\psi(r,z)$ to the hydrostatic
pressure $p(\psi)$ and the the toroidal field $g(\psi)$ (both of which are
assumed to be functions of $\psi$ only), and the currents $I_{k}$ going
through external coils with cross-sectional area $S_{k}$. Posed in free space,
the equation takes the form

$-\nabla\,\cdot\,\left(\frac{1}{\mu
r}\nabla\psi\right)=\left\\{\begin{array}[]{ll}r\frac{d}{d\psi}p(\psi)+\frac{1}{2\,\mu
r}\frac{d}{d\psi}g^{2}(\psi)&\text{ in }\Omega_{p}(\psi)\\\ I_{k}/S_{k}&\text{
in }\Omega_{C_{k}}\\\ 0&\text{ elsewhere. }\end{array}\right.$ (1a) Above,
$\mu$ is the magnetic permeability, $\Omega_{C_{k}}$ denotes the area occupied
by the $k$-th external coil, $\Omega_{p}$ is the plasma confinement region
which is not known a priori and must be determined as a problem unknown,
making this a free boundary problem. A schematic of a cross-section, for
$r>0$, of a tokamak is depicted in Figure 1. The confinement region
$\Omega_{p}$ is characterized as the largest region that contains the magnetic
axis (defined as the point where $\psi$ has a global maximum) and that is
bounded by a closed level set $\psi=constant$. The solution to this free
boundary problem is ubiquitous in nuclear fusion and several computational
codes have been developed over the years (see for instance [14, 20, 24, 25,
28, 31] and references therein).

A common choice in the literature, first proposed in [37], for the free
functions $p(\psi)$ and $g(\psi)$ in the right hand side of (1a) is

$\frac{d}{d\psi}p(\psi)=j_{0}\frac{\beta}{r_{0}}\left(1-\psi_{N}^{\alpha_{1}}\right)^{\alpha_{2}},\qquad\qquad\frac{1}{2}\frac{d}{d\psi}g^{2}(\psi)=j_{0}\mu_{0}r_{0}(1-\beta)\left(1-\psi_{N}^{\alpha_{1}}\right)^{\alpha_{2}},$
(1b)

where $r_{0}$ represents the outer radius of the vacuum chamber and
$\psi_{N}\in[0,1]$ is a normalization of the flux $\psi$ such that
$\psi_{N}=1$ on the plasma boundary $\partial\Omega_{p}$ and $\psi_{N}=0$ at
the magnetic axis. The parameters $\alpha_{1}$ and $\alpha_{2}$ control the
behavior of $\psi$ around the magnetic axis, and $\beta$ measures the ratio
between the hydrostatic pressure and the magnetic pressure in the plasma, and
$j_{0}$ is a normalization factor.

###  2.2 Incorporating uncertainty

In this article, we will consider that the uncertainty in the model (1a) is
concentrated in the values of the currents $I_{i}$ going through the external
coils. As a result, the function $\psi$ (and all quantities derived from it)
are random variables. Obtaining a full description of their probability
density functions might not be possible, but an approximate picture can be
obtained by exploring the parameter space and computing sample approximations
of its expectation and variance. We will model the array of currents as a
$d$-dimensional random variable
$\boldsymbol{\omega}:=(\omega_{1},\ldots,\omega_{d})$, where $d$ is the number
of confinement coils in the reactor, and the $k$-th component of
$\boldsymbol{\omega}$ is the current going through the $k$-th coil. We will
consider that $\boldsymbol{\omega}$ is uniformly distributed around a baseline
vector $\boldsymbol{I}=(I_{1},\ldots,I_{d})$ corresponding to the desired
current values in a deterministic model. We will often refer to
$\boldsymbol{I}$ as either the reference or unperturbed currents. Letting
$\tau>0$ denote the size of the possible perturbation in the current values
(relative to the components of $\boldsymbol{I}$), the vector
$\boldsymbol{\omega}$ is then uniformly distributed over the $d$-dimensional
parameter space

$W:=\prod_{k=1}^{d}\big{[}I_{k}-\tau|I_{k}|,I_{k}+\tau|I_{k}|\big{]}.$ (2)

Since coils are independent of each other, the stochastic random currents
$\left\\{\omega_{k}\right\\}_{k=1}^{d}$ are uncorrelated and the joint density
function of $\boldsymbol{\omega}$ is given by
$\pi\left(\boldsymbol{\omega}\right)=\prod_{k=1}^{d}\pi_{k}\left(\omega_{k}\right)=\prod_{k=1}^{d}\frac{1}{2\tau|I_{k}|}$.
The equilibrium configuration determined by the solution to (1a) is then the
random variable $\psi(x,y,\boldsymbol{\omega})$; we will be primarily
interested in efficiently constructing an approximation to its expected value

$\mathbb{E}\left[\psi(x,y,\boldsymbol{\omega})\right]=\int_{W}\psi(x,y,\boldsymbol{\omega})\pi(\boldsymbol{\omega})d\boldsymbol{\omega},$
(3)

as well as those of some derived quantities such as the plasma boundary, the
location of the x-points, etc.

##  3 Monte Carlo and Multilevel Monte Carlo Finite-Element Methods

We now turn our attention to the numerical approximation of the expected value
(3). Since the location and shape of the plasma boundary depend on the values
of the coil currents, variations of these values could lead to contacts
between the plasma and the wall or even loss of confinement. This fact
translates into a possible non-smoothness of the mapping between coil currents
and the solutions of (1a) which may then cause techniques such as stochastic
collocation to underperform. Moreover, the computational effort associated
with cubature methods scales exponentially with the dimension of the parameter
space, seriously limiting their feasibility for estimating (3). This leads to
the use of Monte Carlo methods, which are agnostic to both the smoothness of
the mapping and the dimensionality of the problem [38], although they have a
slow convergence rate (1/2) that tends to make them costly. This will be
addressed through the use of a multi-level approach.

###  3.1 Monte Carlo Finite-Element method

We will describe the method in terms of a generic solution, $u$, to a PDE
involving stochastic parameters and its finite element approximation $u_{h}$,
where $h$ is the mesh parameter of the discretization. Let
$\\{\boldsymbol{\omega}^{(i)}\\}_{1\leq i\leq N}$ be a set of $N$ realizations
of the random variable $\boldsymbol{\omega}$ giving rise to a sample of $N$
realizations $u^{(i)}=u(\boldsymbol{\omega}^{(i)})$ and
$u_{h}^{(i)}=u_{h}(\boldsymbol{\omega}^{(i)})$ of the exact solution and its
finite element discretization. We will assume that all these functions belong
to a functional space $Z$ endowed with a norm $\|\cdot\|_{Z}$ (see the
Appendix 7.1 for details), and we will consider the standard FEM error
estimate $\|u^{(i)}-u^{(i)}_{h}\|_{Z}\leq C^{(i)}h^{p}$, where $p$ is the
order of the FEM discretization and the constant $C^{(i)}$ depends only on the
problem geometry and the particular values of $\boldsymbol{\omega}^{(i)}$.

The Monte Carlo Finite-Element (MC-FE) estimator $A_{\text{MC}}(u_{h})$ for
$\mathbb{E}(u)$ is defined as the sample mean

$A_{\text{MC}}(u_{h}):=\frac{1}{N}\sum_{i=1}^{N}u_{h}^{(i)}.$ (4)

This estimator is easily shown to be unbiased and to satisfy
$\mathbb{E}(A_{\text{MC}})=\mathbb{E}(u_{h})$. A quantity that serves as a
foundation for examining the spatial and statistical accuracy of the MC-FE
estimator is the mean squared error (MSE) defined as

$\mathcal{E}_{A_{\text{MC}}}^{2}:=\mathbb{E}\left[\left\|\mathbb{E}(u)-A_{\text{MC}}(u_{h})\right\|_{Z}^{2}\right].$

It can be shown (see, for instance [4, Theorem 4.3]) that, for linear
problems, the Monte Carlo estimator accurately approximates the expected value
in the sense that $\mathcal{E}_{A_{\text{MC}}}^{2}\leq K(N^{-1/2}+h^{p})^{2}$,
where the constant $K>0$ depends on the problem geometry and the expected
values of the problem data. The MSE can be decomposed into terms related to
the bias and variance as

$\mathcal{E}_{A_{\text{MC}}}^{2}=\left\|\mathbb{E}(u)-\mathbb{E}(u_{h})\right\|^{2}_{Z}+\mathbb{E}\left[\left\|\mathbb{E}(u_{h})-A_{\text{MC}}(u_{h})\right\|_{Z}^{2}\right]=\left\|\mathbb{E}(u)-\mathbb{E}(u_{h})\right\|^{2}_{Z}+\frac{\mathbb{V}(u_{h})}{N}=\mathcal{E}_{\text{Bias}}^{2}+\mathcal{E}_{\text{Stat}}^{2},$

where $\mathbb{V}(u):=\mathbb{E}[\left\|u-\mathbb{E}(u)\right\|_{Z}^{2}]$ and
$\mathbb{V}(A_{\text{MC}}(u))=\mathbb{V}(u)/N$. The last two terms in the
expression above implicitly define the discretization error
$\mathcal{E}_{\text{Bias}}$ and the sampling (or statistical) error
$\mathcal{E}_{\text{Stat}}$ respectively.

If $\mathbb{E}(u)\neq 0$, the mean squared error can be expressed as a
percentage through normalization by the factor
$\left\|\mathbb{E}(u)\right\|_{Z}^{2}$, leading to the normalized mean squared
error $\widehat{\mathcal{E}}_{A_{\text{MC}}}^{2}$. Since the exact random
variable $u$ is not available, we will approximate the relative mean squared
error (nMSE) by

$\widehat{\mathcal{E}}_{A_{\text{MC}}}^{2}\approx\frac{\left\|\mathbb{E}(u)-\mathbb{E}(u_{h})\right\|^{2}_{Z}}{\left\|\mathbb{E}(u_{h})\right\|_{Z}^{2}}+\frac{\mathbb{V}(u_{h})}{N\left\|\mathbb{E}(u_{h})\right\|_{Z}^{2}}=\widehat{\mathcal{E}}_{\text{Bias}}^{2}+\widehat{\mathcal{E}}_{\text{Stat}}^{2},$
(5)

where $\widehat{\mathcal{E}}_{\text{Bias}}$ and
$\widehat{\mathcal{E}}_{\text{Stat}}$ are relative analogues to the
discretization and statistical errors defined above. If the number of grid
points for the FEM discretization is $M$ then, in two dimensions, it is
standard to assume that
$\widehat{\mathcal{E}}_{\text{Bias}}=\mathcal{O}(M^{-p/2})$. Given a target
tolerance $\epsilon$, the contribution of the statistical and discretization
errors towards the total nMSE can be controlled by requiring that

$\widehat{\mathcal{E}}_{\text{Bias}}^{2}=\mathcal{O}(M^{-p})\leq(1-\theta)\epsilon^{2},\qquad\widehat{\mathcal{E}}_{\text{Stat}}^{2}=\frac{V_{h}}{N}\leq\theta\epsilon^{2},$
(6)

where $\theta\in(0,1)$ is known as the splitting parameter, and
$V_{h}:=\mathbb{V}\left(u_{h}\right)/\left\|\mathbb{E}(u_{h})\right\|_{Z}^{2}$.
This in turn allows us to estimate the sample size $N$ and the number of grid
points $M$ required to achieve the desired tolerance as

$M\geq\left((1-\theta)\epsilon^{2}\right)^{-\frac{1}{p}},\qquad
N=\Biggl{\lceil}\frac{V_{h}}{\theta\epsilon^{2}}\Biggr{\rceil}.$ (7)

Assuming the average cost to obtain one sample (i.e. to solve (1a) for one
value of the coil-currents) is $C=\mathcal{O}(M^{c}$) for some $c>0$, the
total computational cost of the MC-FE estimator can be estimated as

$C(A_{\text{MC}})=\mathcal{O}(NM^{c})=\mathcal{O}\left(\epsilon^{-2-\frac{2c}{p}}\right).$
(8)

###  3.2 Multilevel Monte Carlo Finite-Element method

MLMC reduces the computational cost associated with sampling—which in our case
involves the numerical solution of a non-linear PDE in a target computational
mesh—by approximating the expectation of the quantity of interest on the
finest mesh by a sequence of control variates on a set of coarse grids that
are cheaper to compute [13]. Using the linearity of expectation, the MLMC
estimator expresses the quantity of interest on the finest spatial grid,
$\mathbb{E}\left(u_{h}\right)$, by a telescoping sum involving the numerical
approximations of $u$ on coarser grids. Consequently, MLMC’s workload is
shifted from the fine grid to coarser grids, making it more efficient than MC
[8]. Cf. e.g. [35] for alternative ways to reduce the costs of MC methods.

To construct meshes that are easy to describe for both uniform and adaptive
mesh refinement, we will characterize them using the number of grid points
rather than the mesh size. We will refer to $\ell=0,\dots,L$ as the level of a
mesh $\mathcal{T}_{\ell}$ containing $\\{M_{\ell}\\}$ grid points. We will
then consider a sequence of meshes $\mathcal{T}_{0},\ldots,\mathcal{T}_{L}$
with increasing resolution so that $\\{M_{\ell}\\}_{0\leq\ell\leq L}$ defines
an increasing sequence and $\mathcal{T}_{L}$ is the finest mesh, and we will
denote by $u_{\ell}$ the approximation of $u$ on the mesh
$\mathcal{T}_{\ell}$.

The expectation of the function $u$ can be approximated by the expectation of
the finest approximation $u_{L}$. Since
$u_{L}=u_{0}+(u_{1}-u_{0})+(u_{2}-u_{1})+\cdots+(u_{L}-u_{L-1})$, this
approximation is given by

$\mathbb{E}(u)\approx\mathbb{E}(u_{L})=\mathbb{E}(u_{0})+\sum_{\ell=1}^{L}\mathbb{E}(u_{\ell}-u_{\ell-1})=\sum_{\ell=0}^{L}\mathbb{E}(Y_{\ell}),$
(9)

a telescoping sum, where each of the terms

$Y_{0}:=u_{0}\quad\text{ and }\quad
Y_{\ell}:=u_{\ell}-u_{\ell-1}\quad(\text{for }\ell\geq 1)$ (10)

can be regarded as a correction of the coarsest approximation $u_{0}$. If each
term $\mathbb{E}(Y_{\ell})$ is estimated by gathering $N_{\ell}$ samples at
level $\ell$ and computing the sample expectations

$\mathbb{E}(Y_{0})\approx\widehat{Y}_{0}:=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}u_{0}^{(i)},\qquad\mathbb{E}(Y_{\ell})\approx\widehat{Y}_{\ell}:=\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}\left(u_{\ell}^{(i)}-u_{\ell-1}^{(i)}\right)\quad(\text{
for }\ell\geq 1),$

then the MLMC-FE estimator at level $L$ will be unbiased and can be written as

$A_{\text{MLMC}}(u_{L}):=\sum_{\ell=0}^{L}\widehat{Y}_{\ell}=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}u_{0}^{(i)}+\sum_{\ell=1}^{L}\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}\left(u_{\ell}^{(i)}-u_{\ell-1}^{(i)}\right).$
(11)

Recalling that $\mathbb{E}(\widehat{Y}_{\ell})=\mathbb{E}(Y_{\ell})$ and
$\mathbb{V}(\widehat{Y}_{\ell})=\mathbb{V}(Y_{\ell})/N_{\ell}$ we conclude
that, for the MLMC-FE estimator, it follows that
$\mathbb{E}\left(A_{\text{MLMC}}\right)=\sum_{\ell=0}^{L}\mathbb{E}(\widehat{Y}_{\ell})=\mathbb{E}(u_{L})$
and
$\mathbb{V}\left(A_{\text{MLMC}}\right)=\sum_{\ell=0}^{L}\mathbb{V}(\widehat{Y}_{\ell})=\sum_{\ell=0}^{L}\mathbb{V}(Y_{\ell})/N_{\ell}$.

As for MC-FE, the mean squared error $\mathcal{E}_{A_{\text{MLMC}}}^{2}$ can
be split into contributions from bias and variance as

$\mathcal{E}_{A_{\text{MLMC}}}^{2}=\mathbb{E}\left[\left\|\mathbb{E}(u)-A_{\text{MLMC}}(u_{L})\right\|_{Z}^{2}\right]=\left\|\mathbb{E}(u)-\mathbb{E}(u_{L})\right\|^{2}_{Z}+\sum_{\ell=0}^{L}\frac{\mathbb{V}\left(Y_{\ell}\right)}{N_{\ell}}=\mathcal{E}_{\text{Bias}}^{2}+\mathcal{E}_{\text{Stat}}^{2}.$

Similarly, using $\mathbb{E}(u)\approx\mathbb{E}(u_{L})$, the normalized mean
squared error $\widehat{\mathcal{E}}_{A_{\text{MLMC}}}^{2}$ can be
approximated by

$\widehat{\mathcal{E}}_{A_{\text{MLMC}}}^{2}\approx\frac{\left\|\mathbb{E}(u)-\mathbb{E}(u_{L})\right\|^{2}_{Z}}{\left\|\mathbb{E}(u_{L})\right\|_{Z}^{2}}+\sum_{\ell=0}^{L}\frac{V_{\ell}}{N_{\ell}}=\widehat{\mathcal{E}}_{\text{Bias}}^{2}+\widehat{\mathcal{E}}_{\text{Stat}}^{2},$
(12)

where

$V_{\ell}:=\mathbb{V}\left(Y_{\ell}\right)/\left\|\mathbb{E}(u_{L})\right\|_{Z}^{2}.$
(13)

As above, the parameter $\theta\in(0,1)$ can be used to split the
contributions of the two components of the error by requiring
$\widehat{\mathcal{E}}_{\text{Bias}}^{2}\leq(1-\theta)\epsilon^{2}$ and
$\widehat{\mathcal{E}}_{\text{Stat}}^{2}\leq\theta\epsilon^{2}$ where
$\epsilon$ is a predetermined tolerance such that
$\widehat{\mathcal{E}}_{A_{\text{MLMC}}}^{2}\leq\epsilon^{2}$.

Note that each of the terms $Y_{\ell}^{(i)}:=u_{\ell}^{(i)}-u_{\ell-1}^{(i)}$
appearing in (11) requires the approximation of $u^{(i)}$ on adjacent
refinement levels using the same value of the parameter
$\boldsymbol{\omega}^{(i)}$. However, for FEM discretization, the numerical
implementation of this term does not require the solution of the PDE on the
two grids $\mathcal{T}_{\ell}$ and $\mathcal{T}_{\ell-1}$. The realization on
the coarse grid can be obtained from the one on the fine grid by either
Galerkin projection or interpolation. Projection is the more accurate choice,
but it requires the solution of a system involving the mass matrix for each
realization. Thus we prefer interpolation, for which costs are minimal and
which is of second-order accuracy in our application. To avoid introducing
correlation across discretization levels, none of the samples involved in the
computation of $Y_{\ell}$ is reused for the finer level $\ell+1$. That is,
sampling is done such that, for $n\neq m$, the estimates $Y_{n}$ and $Y_{m}$
are uncorrelated. However, for any particular $Y_{\ell}$ the strong
correlation between $u_{\ell}^{(i)}$ and $u_{\ell-1}^{(i)}$ makes the variance
of the correction terms much smaller than the variance of the approximation
$u_{L}$ in the finest mesh, further improving the statistical approximation.

To quantify the computational effort of the MLMC-FE estimator, let
$M_{\ell,i}$ be the number of grid points for the $i$-th sample on mesh level
$\ell$. We will assume that the computational cost to obtain one sample of
$u_{\ell}^{(i)}$ is
$C_{\ell,i}:=C(u_{\ell}^{(i)})=\mathcal{O}(M_{\ell,i}^{c})$, where the
exponent $c>0$ depends on the solver, and will denote the cost of computing
the correction term $Y_{\ell}^{(i)}$ by
$C_{\ell,i}:=C(Y_{\ell}^{(i)})=\mathcal{O}(M_{\ell,i}^{c})$ for $\ell\geq 0$
and $M_{-1,i}=0$. For a nonlinear problem like the one at hand, the particular
realization $\boldsymbol{\omega}^{(i)}$ will influence the cost. We will
consider the average cost to be of the form

$C_{\ell}=\mathcal{O}(M_{\ell}^{c})$ (14)

and use this to estimate the total cost as
$C(A_{\text{MLMC}})=\sum_{\ell=0}^{L}N_{\ell}C_{\ell}=\sum_{\ell=0}^{L}N_{\ell}\cdot\mathcal{O}\left(M_{\ell}^{c}\right).$
Using the method of Lagrange multipliers, it is shown in [19] that this total
cost, a function of $N_{\ell}$, can be minimized subject to the constraint
$\widehat{\mathcal{E}}_{\text{Stat}}^{2}\leq\theta\epsilon^{2}$ by the choice

$N_{\ell}=\Biggl{\lceil}\frac{1}{\theta\epsilon^{2}}\sqrt{\frac{V_{\ell}}{C_{\ell}}}\sum_{k=0}^{L}\sqrt{V_{k}C_{k}}\Biggr{\rceil}.$
(15)

With this expression for $N_{\ell}$, the optimal total cost for the MLMC-FE
estimator is

$C(A_{\text{MLMC}})\leq\frac{1}{\theta\epsilon^{2}}\left(\sum_{\ell=0}^{L}\sqrt{V_{\ell}C_{\ell}}\right)^{2}+\sum_{\ell=0}^{L}C_{\ell}.$
(16)

This analysis of computational cost is established in full generality in [19,
Theorem 1].

The formula (15) suggests an iterative procedure for the approximation of
$\mathbb{E}(u)$. Starting from a computational mesh $\mathcal{T}_{0}$, gather
an initial number $\widetilde{M}_{0}$ of samples $u_{0}^{(i)}$ and estimate
$\widehat{\mathcal{E}}_{\text{Bias}}$, $\widehat{\mathcal{E}}_{\text{Stat}}$,
and $V_{0}$. If $\widehat{\mathcal{E}}_{\text{Stat}}$ is larger than the
prescribed tolerance, use (15) to update $\widetilde{M}_{0}$ and gather
additional samples; if $\widehat{\mathcal{E}}_{\text{Bias}}$ is above the
prescribed tolerance, then add an additional level of spatial refinement. The
process continues adding discretization levels and collecting additional
and/or samples until both $\widehat{\mathcal{E}}_{\text{Bias}}$ and
$\widehat{\mathcal{E}}_{\text{Stat}}$ fall below the required tolerance, at
which point $\mathbb{E}(u)$ is approximated using equation (9). This simple
algorithm, however, presents one challenge: the term $V_{\ell}$ in equation
(15) requires the computation of the term

$\mathbb{V}(Y_{\ell})=\frac{1}{N_{\ell}-1}\left(\sum_{i=1}^{N_{\ell}}\left\|Y_{\ell}^{(i)}\right\|_{Z}^{2}-\frac{1}{N_{\ell}}\left\|\sum_{i=1}^{N_{\ell}}Y_{\ell}^{(i)}\right\|_{Z}^{2}\right).$
(17)

However, the estimated sample sizes $\\{N_{\ell}\\}$ are available only for
pre-existing discretization levels, hence whenever an additional mesh
refinement is needed, $N_{L+1}$ cannot be approximated by (15) as this formula
uses $V_{L+1}$ to compute $N_{L+1}$. This inconvenience can be overcome by
noting that

$\mathbb{V}(u-u_{\ell})=\mathbb{E}\left[\left\|u-u_{\ell}\right\|_{Z}^{2}\right]-\left\|\mathbb{E}\left(u-u_{\ell}\right)\right\|_{Z}^{2}\leq\mathbb{E}\left[\left\|u-u_{\ell}\right\|_{Z}^{2}\right].$
(18)

Hence, the variance $\mathbb{V}(u-u_{\ell})$ can be estimated by the
expectation of the squared discretization error
$\left\|u-u_{\ell}\right\|_{Z}^{2}$. For a uniformly refined grid, we can
resort to a standard a priori error estimate and assume that
$\mathbb{E}[\left\|u-u_{\ell}\right\|_{Z}^{2}]=\mathcal{O}(M_{\ell}^{-b_{1}})$
for some $b_{1}>0$. It then follows that
$V_{\ell}=\mathcal{O}(M_{\ell}^{-b_{1}})$ and following [40] we can then
approximate $V_{L+1}$ in terms of the known variance $V_{L}$ by

$V_{L+1}=\left(M_{L+1}/M_{L}\right)^{-b_{1}}V_{L}.$ (19)

###  3.3 Adaptive multilevel Monte Carlo Finite-Element method

In view of the benefits of approximating the quantity of interest across a
sequence of increasingly finer meshes, and with the goal of further reducing
the computational cost associated with reducing the bias associated with the
numerical discretization, it is natural to focus the refinement only on those
parts of the mesh where the error is concentrated. Our goal is then to,
starting from a computational mesh $\mathcal{T}_{0}$, generate a family of
adaptively refined meshes $\left\\{\mathcal{T}_{\ell}\right\\}_{0\leq\ell\leq
L}$ that will produce better approximations of $\mathbb{E}(u)$ than the ones
resulting from consecutive uniform refinements of the initial grid. With this
goal in mind, the use of an a posteriori error estimator to guide the
construction of the family of meshes has been proposed in the context of
multilevel Monte Carlo methods [11, 26, 27, 32].

Residual error estimator. A key ingredient in an adaptive solver is a local
error estimator. In our case, for each element $K$ on the mesh
$\mathcal{T}_{\ell}$, we will use the simple residual-based a posteriori error
indicator

$\eta_{K,\ell}(\boldsymbol{\omega}):=h_{K}^{2}\left\|\,\nabla\cdot\left(\tfrac{1}{\mu
r}\nabla
u_{\ell}(\boldsymbol{\omega})\right)-f(u_{\ell}(\boldsymbol{\omega}))\,\right\|_{K}+h_{K}^{3/2}\left\|\,\left[\\!\left[\tfrac{1}{\mu
r}\nabla
u_{\ell}(\boldsymbol{\omega})\cdot\boldsymbol{n}\right]\\!\right]\,\right\|_{\partial
K\setminus\partial\Omega},$ (20)

where $\partial K$ is the boundary of the element $K$, $h_{K}$ is the diameter
of $K$, $\boldsymbol{n}$ is the outward unit normal to the element $K$,
$[\\![\cdot]\\!]$ denotes the jump across the edge of an interior element, and
$f$ is the source term defined piecewise on the right-hand side of (1a).
Following [11, 32], we will further define the mean local and mean global
error estimators respectively as

$\eta_{K,\ell}:=\mathbb{E}\left(\eta_{K,\ell}(\boldsymbol{\omega})\right)\quad\text{
and }\quad\eta_{\ell}^{2}:=\sum_{K\in\mathcal{T}_{\ell}}\eta_{K,\ell}^{2}.$
(21)

For linear deterministic problems, estimators of this form can be shown to be
such that there are constants $C_{1},C_{2}>0$ such that
$C_{1}\eta_{\ell}\leq\|u-u_{\ell}\|_{Z}\leq C_{2}\eta_{\ell}$. Therefore, the
error estimator will accurately locate the regions of high error density and
will decay at the same rate as the true error [34]. The global error can then
be approximated by adding the local estimators over the entire triangulation.
Using these error estimators, the adaptive analogue of (18) can be written as

$\mathbb{V}(u-u_{\ell})\leq\mathbb{E}\left[\left\|u-u_{\ell}\right\|_{Z}^{2}\right]\approx\eta_{\ell}^{2},$
(22)

which then leads to the following adaptive analogue of the extrapolation
formula (19)

$V_{L+1}=\left(\eta_{L+1}/\eta_{L}\right)^{2}V_{L}.$ (23)

This estimate can then be used in combination with (15) to obtain an update
for the sample size required at each adaptive level.

Adaptive solution cycle. With these definitions in place, we can then describe
our strategy, which follows the “SOLVE $\rightarrow$ ESTIMATE $\rightarrow$
MARK $\rightarrow$ REFINE” paradigm familiar from deterministic adaptive
solvers [7], as:

  1. 1.

Solve: Starting from a fixed number of samples, the problem (1a) is solved on
the initial mesh $\mathcal{T}_{0}$.

  2. 2.

Estimate: The local mean error estimator is approximated from the sample
gathered.

  3. 3.

Mark: The set $\mathcal{M}_{\ell}$ containing the smallest possible number of
elements in $\mathcal{T}_{0}$ satisfying

$\sum_{K\in\mathcal{M}_{\ell}}\eta_{K,\ell}^{2}\geq\zeta\eta_{\ell}^{2},$ (24)

for some predetermined value $\zeta\in[0,1]$ is marked for refinement—this
marking strategy is known as  Dörfler marking in the adaptive finite element
community [10].

  4. 4.

Refine: The elements marked are refined in such a way that the resulting
triangulation $\mathcal{T}_{\ell+1}$ is shape-regular and conforming. Efforts
should be made to make sure that the growth of the number of elements is kept
at bay. In our case, we used the implementation given in [16] of the
algorithms described in [6, 15, 39].

The steps above are repeated until the error estimator $\eta$ falls below a
certain predetermined value.

A notion of mesh level. For uniformly refined grids, the notion of mesh level
is natural: starting from a mesh $\mathcal{T}_{\ell}$, one step of uniform
refinement decreases the mesh parameter $h$ across the grid by a factor of
$1/2$; the resulting mesh is said to have level $\ell+1$ and is denoted by
$\mathcal{T}_{\ell+1}$. For adaptively refined grids, where the mesh parameter
is not constant through the grid, the notion of the level does not come so
naturally. We will use the fact that, for a uniform refinement, the numerical
error decays by a factor of $(1/2)^{p}$, (where $p$ is the order of the FEM
solver) with each successive level to extend the notion of mesh level to
adaptively refined grids.

Consider a numerical approximation $u_{\ell}$ obtained on a mesh
$\mathcal{T}_{\ell}$ with an associated error estimation given by
$\eta_{\ell}$. We will say that a mesh has level $\ell+1$ and will denote it
by $\mathcal{T}_{\ell+1}$ if it was obtained from $\mathcal{T}_{\ell}$ by
cycling over the adaptive loop using the value $(1/2)^{p}\eta_{\ell}$ as the
stopping tolerance. In other words, we will say that an adaptively refined
mesh has level $\ell+1$ if it produces a numerical solution with an error
$(1/2)^{p}$ times smaller than one with level $\ell$, just like in the uniform
case. We will refer to $q:=(1/2)^{p}$ as the decay factor. In terms of
discretization accuracy, after $\ell$ steps of adaptive refinement, an
adaptively refined mesh with level $\ell$ will have an associated error
estimation $\eta_{\ell}=q^{\ell}\eta_{0}$, where $\eta_{0}$ corresponds to the
error estimation at the initial mesh. In our numerical experiments, since the
convergence rate of the piecewise linear solver is 2 (when measured in the
$L^{2}$ norm), we shall use a decay factor $q=1/4$ to define our adaptively
refined meshes.

Deterministic adaptive grids. Ideally, in the stochastic setting, all the
error estimations collected from the totality of samples would be used to
drive the adaptive refinement forward and build an optimal set of meshes at
every level. However, due to the iterative nature of the algorithm arising
from (15), the optimal mesh at every level would have to be corrected with
every new batch of samples and the solutions corresponding to all previous
realizations $\boldsymbol{\omega}^{(i)}$ would have to be recomputed. The
computational cost of re-sampling in this manner quickly becomes impractical.

Instead, following [26, 27, 33], we will construct a sequence of deterministic
adaptive grids with partial knowledge about $\mathbb{E}(u)$ as follows.
Starting from a sample $\\{\boldsymbol{\omega}^{(i)}\\}_{1\leq i\leq N}$
(where $N$ is small and arbitrarily chosen) and a mesh $\mathcal{T}_{0}$, the
PDE is solved and the local error is estimated for every solution
$u_{0}^{(i)}$, resulting in $N$ local error estimators
$\\{\eta_{K,0}(\boldsymbol{\omega}^{(i)})\\}_{1\leq i\leq N}$. The mean local
and global error estimators $\eta_{K,0}$ and $\eta_{0}$ are then approximated
by the sample means of the individual estimators. Using this approximation of
$\eta_{K,0}$, the mesh $\mathcal{T}_{0}$ is refined. This process is continued
until the approximated mean error estimator satisfies $\eta\leq q\eta_{0}$;
the resulting mesh is stored and labeled as $\mathcal{T}_{1}$ (mesh level 1).
The previous steps are repeated until a target number of meshes
$\\{\mathcal{T}_{\ell}\\}_{0\leq\ell\leq L}$ have been generated. The process
is described in Algorithm 1. Since the family of meshes produced is
constructed using random samples of $\boldsymbol{\omega}$, they approximately
reduce the error for the approximate expectation $\mathbb{E}(u_{h})$ by a
factor of $q$ with every increasing level. This family of meshes is then kept
fixed during the MLMC run.

1

Input: Initial mesh $\mathcal{T}_{0}$, maximum mesh level $L$, element marking
percentage $\zeta\in(0,1)$, sample size $N$, error decay factor $q\in(0,1)$.

Output: Adaptive mesh set $\\{\mathcal{T}_{\ell}\\}_{\ell=0}^{L}$.

2

3for  _$\ell=1,\cdots,L$_ do

4 Refine = TRUE.

5 Set $j=1$ and $\mathcal{T}^{(j)}=\mathcal{T}_{\ell-1}$.

6 while  _Refine_ do

7

8 for  _$i=1,\cdots,N$_ do

9 (i) Draw a random sample $\boldsymbol{\omega}^{(i)}$ from $W$.

10 (ii) Obtain $u^{(i)}$ by solving the free boundary problem on on
$\mathcal{T}^{(j)}$ using $\boldsymbol{\omega}^{(i)}$.

11 (iii) Calculate and accumulate the local and global error estimators
$\eta_{K,j}(\boldsymbol{\omega}^{(i)})$ and
$\eta_{j}(\boldsymbol{\omega}^{(i)})$.

12 Compute the mean estimators $\eta_{K,j}$ and $\eta_{j}$ from the
accumulated samples.

13 Mark & refine the mesh $\mathcal{T}^{(j)}$ according to $\eta_{K,j}$ to
obtain $\mathcal{T}^{(j+1)}$.

14 if  _ $\eta_{j}\leq q\eta_{1}$_ then

15 Store $\mathcal{T}^{(j)}$ as $\mathcal{T}_{\ell}$.

16 Refine = FALSE.

17

18 $j=j+1$

19

Algorithm 1 Generate adaptive deterministic mesh set

1

Input: Initial mesh level $L=0,$ sequence of mesh available
$\\{\mathcal{T}_{\ell}\\}_{\ell\geq 0}$, root nMSE $\epsilon$,
$\zeta\in(0,1)$, initial sample size
$N_{\text{old}}=\left\\{N_{\ell}\right\\}_{\ell=0}^{L}$, counter $j=1$, sample
size corrections $dN_{\ell}^{j}=\\{N_{\ell}\\}_{\ell=0}^{L}$.

Output: $\left\\{N_{\ell}\right\\}_{\ell=0}^{L}$,
$\\{V_{\ell}\\}_{\ell=0}^{L}$, $A_{\text{MLMC}}$.

2 while  _$\sum_{\ell}(dN_{\ell}^{j}) >0$_ do

3 for  _$0\leq\ell\leq L$_ do

4

5 for  _$i=1,\cdots,dN_{\ell}^{j}$_ do

6 Solve the free boundary problem on $\mathcal{T}_{\ell}$ to get
$u_{\ell}^{(i)}$ for the $i$-th sample.

7

8 Approximate $\\{V_{\ell}\\}_{\ell=0}^{L}$ by (17).

9 Update the sample size estimation $\left\\{N_{\ell}\right\\}_{\ell=0}^{L}$
by (15).

10 $j=j+1$.

11 $dN_{\ell}^{j}=\left\\{N_{\ell}\right\\}_{\ell=0}^{L}-N_{\text{old}}.$

12 $N_{\text{old}}=\left\\{N_{\ell}\right\\}_{\ell=0}^{L}$.

13 if  _$\sum_{\ell}(dN_{\ell}^{j})=0$_ then

14

15 if  _The a posteriori error estimator falls below
$\sqrt{1-\theta}\epsilon\left\|\mathbb{E}(u)\right\|_{Z}$, _ then

16 Compute $A_{\text{MLMC}}$ by (11) and terminate the loop.

17 else

18 $L=L+1$.

19 Approximate $V_{L}$ by (19) and compute
$\left\\{N_{\ell}\right\\}_{\ell=0}^{L}$ and go to Step 1.

20

21

Algorithm 2 Multilevel Monte Carlo Finite-Element method

##  4 Numerical experiments

We now demonstrate and evaluate the performance of the methods presented for
the Grad-Shafranov free boundary problem. We examine the efficiency of the
three simulation approaches, comparing the CPU times (or computational cost)
and the accuracy of some geometric descriptors that are generated from the
approximation of $\mathbb{E}(\psi)$ obtained with each of the techniques.
Following [14], we consider an ITER geometry with 12 coils and a “baseline”
vector of target current intensities $\boldsymbol{I}$ given by

${\begin{array}[]{llll}I_{1}=-1.4\times 10^{6}A,&I_{2}=-9.5\times
10^{6}A,&I_{3}=-2.0388\times 10^{7}A,&I_{4}=-2.0388\times 10^{7}A,\\\
I_{5}=-9\times 10^{6}A,&I_{6}=3.564\times 10^{6}A,&I_{7}=5.469\times
10^{6}A,&I_{8}=-2.266\times 10^{6}A,\\\ I_{9}=-6.426\times
10^{6}A,&I_{10}=-4.82\times 10^{6}A,&I_{11}=-7.504\times
10^{6}A,&I_{12}=1.724\times 10^{7}A.\end{array}}$ (25)

We will refer to these values as the reference currents. The profiles for $p$
and $g$ on the right hand side of (1a) follow the form given in (1b) with the
specific values $r_{0}=6.2m,\,\beta=0.5978,\,\alpha_{1}=2$, and
$\alpha_{2}=1.395$, and $j_{0}=1.3655\times 10^{6}A/m^{2}$. The reactor and
coil array geometries follow the ones described in [2]. In our experiments, we
will take the vector of current intensities to be subject to uncertainty
modeled as a uniformly distributed perturbation of magnitude $\tau=2\%$
centered around the reference values above.

###  4.1 Experiment description

In this section, we present numerical results comparing the three approaches –
MC-FE, uniform MLMC-FE, and adaptive MLMC-FE. For the solution of (1), we used
the finite element-based solver FEEQS.m [23] developed by Holger Heumann and
collaborators as a lightweight Matlab implementation of the code CEDRES++ [14,
24]. The code implements a piecewise linear finite element discretization of a
weak formulation of (1) and employs a globalized variation of Newton’s method
to resolve the nonlinearity. (The stopping threshold for the relative residual
was set to $5\times 10^{-11}$.) For the solution of the perturbed problems,
the initial iterate of Newton’s method was taken to be the solution
corresponding to the reference currents $\boldsymbol{I}$. All tests used the
splitting parameter $\theta=0.5$ in (6). The user-specified tolerances for the
normalized mean squared error range from $\epsilon=2\times 10^{-4}$ to
$8\times 10^{-3}$. Experiments were conducted using MATLAB R2022a on a System
76 Thelio Major with 256GB RAM and a 64-Core @4.6 GHz AMD Threadripper 3
Processor.

To produce an estimate of the number of samples required on each
discretization level, equation (15) requires the knowledge of two problem-
dependent parameters: the power $c$ appearing in the estimate for the
computational cost (14), and the normalized variance of the correction terms
$V_{\ell}$ as defined in (13). The normalization factor
$\|\mathbb{E}(u_{L})\|_{Z}$ in (13) was estimated on the finest uniform mesh
level ($\ell=5$) to be approximately $8.57\times 10^{-1}$. To estimate the
value of $c$, 100 random currents are sampled for different mesh sizes
$M_{\ell}$, the processing times required to obtain the solutions are averaged
for each mesh size, and $c$ is estimated through a regression. Figure 2 shows
the behavior of the average cost as a function of the mesh size $M_{\ell}$;
from the data displayed, the power law is estimated to be $c\approx 1.09$.
Note that this cost estimate is based on Matlab timings and not on the
complexity analysis of standard linear solution algorithms. The same samples
are also used to estimate the sample means, $\mathbb{E}(Y_{\ell})$ or
$\mathbb{E}(u_{h})$, and variances, $\mathbb{V}(Y_{\ell})$ or
$\mathbb{V}(u_{h})$ dynamically using Welford’s algorithm [45]. As the new
samples are gathered, the mean $m_{w}$ and proxy for the variance, $s_{w}$,
are updated using the following formulas for the $i$-th sample, with
$m_{w}^{(0)}=0,s_{w}^{(0)}=0$:

$m_{w}^{(i)}=m_{w}^{(i-1)}+\frac{u^{(i)}-m_{w}^{(i-1)}}{i},\quad
s_{w}^{(i)}=s_{w}^{(i-1)}+\left\langle
u^{(i)}-m_{w}^{(i-1)},u^{(i)}-m_{w}^{(i)}\right\rangle.$

Using sample size $i$, the variance is then given by
$V^{(i)}=s_{w}^{(i)}/(i-1)$.

Level $\ell$ | 0 | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
Uniform | 2685 | 8019 | 30449 | 120697 | 484080 | 1934365
Adaptive | 2685 | 6090 | 25099 | 103968 | 411913 | 1552282

Table 1: The number of grid points $M_{\ell}$ for both geometry-conforming
uniform and adaptive ($q=1/4$) meshes as the mesh levels increase from 0 to 5.

To perform MLMC-FE simulations, the user typically defines and generates a
sequence of spatial grids, where, given a tolerance $\epsilon$, the fineness
of the grid is determined by the requirement that the discretization error
($\widehat{\mathcal{E}}_{\text{Bias}}^{2}$ in (6)) or an estimate of it be
small enough. In this study, we generated two types of grids, a set of
geometry-conforming uniformly refined grids, and a set of adaptively refined
grids constructed using the strategy presented in Section 3.3.111The domain
components contain curved boundaries, which we handled by treating them as
polygonal structures. The mesh generation entails identifying the curved
boundaries using piecewise splines and interpolating along these splines with
grids of varying fineness. For the uniformly refined grid, we generated a
total of six levels of grids. We note that, due to the increasing accuracy of
the spline approximation of the curved boundaries, these meshes are not a
direct refinement of each other. Instead, each level is characterized by a
mesh parameter $h_{\ell}$ being roughly half of the preceding mesh as
determined by the Triangle mesh generator [43]. The adaptive refinement
strategy began with the coarsest mesh from the uniform family and applied the
weighted $L_{2}$ a posteriori error estimator specified in (20) and $q=1/4$ to
reflect a similar error decay as for uniform refinement, also using Triangle
to generate the desired adaptively refined meshes. The number of grid points
for each of these methods on different grid levels is shown in Table 1. Note
that the grid sizes for the adaptive meshes are not dramatically smaller than
for uniform meshes, which suggests that the solution, as a function of the
spatial variables, is fairly smooth.

---
Figure 2: Mean CPU time to compute 100 realizations of $u_{\ell}$, as a
function of the number of grid points $M_{\ell}$, plotted on a logarithmic
scale. The fitted curve indicates that the computational cost $C_{\ell}$
behaves approximately like $M_{\ell}^{1.09}$.

###  4.2 Computational cost

Figure 3 shows a variety of computational results, including the error
estimator, $V_{\ell}$, and CPU time for the two versions of MLMC (and times
for full MC). To investigate the convergence behavior of the discretization
error, we calculate the a posteriori error estimator for both uniform and
adaptive meshes in the same experiment to obtain an estimate of $c$ for
$C_{\ell}$ before conducting the simulations. The results are displayed in the
top left plot in Figure 3, with a dashed line showing a least square fit,
indicating that the discretization error of both methods exhibits an
asymptotic rate of $\mathcal{O}(M^{-1})$ (or $p\approx 2$). The similar
convergence rate further indicates that the solution to the problem is smooth
and the error is equidistributed, rendering the adaptive strategy comparable
to uniform mesh refinement. Note that the estimated error is used for variance
extrapolation in (23) during the MLMC simulations.

The top right plot of Figure 3 shows the behavior of $V_{\ell}$ of (13) for
both uniform and adaptive MLMC-FE methods with $\epsilon=2\times 10^{-4}$,
using six levels of meshes. It can be seen that both methods demonstrate a
decreasing trend in the values of $V_{\ell}$ as the mesh resolution increases,
with a power law decay characterized by $b\approx 2$ in the least square fit.
But there is a regime for a small number of grid points where the asymptotic
behavior of the adaptive method is not evident, in contrast to the behavior of
the uniform method. As the meshes get finer, the plots of $V_{\ell}$ for the
two methods are close to being parallel.

The computational effort for uniform MLMC-FE and MC-FE scales as
$\mathcal{O}(\epsilon^{-2})$ and $\mathcal{O}(\epsilon^{-3})$ respectively, as
indicated by the slopes of the least square fitting lines for the red and
yellow curves. This observation is consistent with the theoretical cost
predictions in Theorem 1 (with $b>c$) and (8). Theorem 1 also indicates that
the majority of computational work is performed on coarse grids. Table 2 shows
the sample sizes obtained from (7) for MC-FE and (15) for MLMC-FE, further
demonstrating a decrease in $N_{\ell}C_{\ell}$ as $\ell$ increases for the
multilevel methods. We also found that the computational cost associated with
the smallest tolerance $\epsilon=2\times 10^{-4}$ is so large that we were
unable to generate MC-FE results on a fine mesh ($\ell=5$) with a large sample
size. In contrast, both versions of MLMC-FE successfully generated results
with this tolerance.222Although we could not directly generate results for MC-
FE for $\epsilon=2\times 10^{-4}$, we could estimate the costs. In particular,
we found that the variance $V_{h}$ is close to constant across mesh levels.
Consequently, we used (7) to estimate the number of required samples as 8000
in Table 2, approximately four times the number required for $\epsilon=4\times
10^{-4}$. This number was multiplied by the mean CPU time observed for the
computations for Figure 2 (120.3 seconds, the largest entry appearing in the
figure) to give the estimated total CPU time in Table 3. Table 3 gives
quantitative data on the costs in CPU time for the three methods, as well as
the speedups achieved by the two multilevel methods. It can be seen that for
small tolerance $\epsilon$, both these methods reduce the CPU times
dramatically, with many examples of speedups greater than a factor of 60 and a
best-case speedup of approximately 200.

As seen in the bottom left plot of Figure 3, the uniform and adaptive versions
of MLMC-FE have similar computational costs of $\mathcal{O}(\epsilon^{-2})$,
as evidenced by the similar decay rate of the error estimator and (22).
According to (16), the slightly smaller magnitude of the error estimator for
the adaptive MLMC-FE suggests a smaller (or comparable) computational cost in
the asymptotic regime. However, when $\epsilon=4\times 10^{-4}$, the adaptive
MLMC-FE method requires approximately twice as much CPU time ($1.79\times
10^{3}$ seconds) compared to the uniform MLMC-FE approach ($9.29\times 10^{2}$
seconds) due to a notable increase in $V_{\ell}$ around $M_{\ell}=10^{4}$.
This also causes the speedups achieved using adaptive refinement to be
somewhat smaller than for uniform refinement. Thus, the traditional advantage
of adaptive mesh refinement is not clearly present. We also attribute this to
the apparent smoothness of the solution. We will demonstrate some advantages
of the adaptive strategy in Section 4.3.

The bottom right plot of Figure 3 shows that the nMSE tolerance $\epsilon$ of
MC-FE approach declines at $\mathcal{O}(N^{-0.51})$, which is consistent with
the well-known square root convergence rate. This rate holds since
$\mathbb{V}(u_{h})$ remains nearly constant among all levels.

|
---|---
|
Figure 3: Top left: weighted $L_{2}$ error (with weight $\mu x$) of estimator
$C\eta_{\ell}$ vs. number of grid points $M_{\ell}$ plot. Top right:
normalized variance $V_{\ell}$ vs. number of grid points $M_{\ell}$ plot.
Bottom left: CPU time in seconds vs. tolerance $\epsilon$. Bottom right: Monte
Carlo convergence rate estimate with tolerance $\epsilon$ vs. sample size $N$.
This plot is generated from Table 2.

| Level $\ell$
---|---
$\epsilon$ | 0 | 1 | 2 | 3 | 4 | 5
$8\times 10^{-3}$ |  |  | 5 |  |  |
$6\times 10^{-3}$ |  |  | 7 |  |  |
$4\times 10^{-3}$ |  |  |  | 22 |  |
$2\times 10^{-3}$ |  |  |  | 83 |  |
$10^{-3}$ |  |  |  |  | 322 |
$8\times 10^{-4}$ |  |  |  |  | 527 |
$6\times 10^{-4}$ |  |  |  |  | 869 |
$4\times 10^{-4}$ |  |  |  |  | 1980 |
$2\times 10^{-4}$ |  |  |  |  |  | 8000∗
| Level $\ell$
---|---
$\epsilon$ | 0 | 1 | 2 | 3 | 4 | 5
$8\times 10^{-3}$ | 10 | 2 | 2 |  |  |
$6\times 10^{-3}$ | 12 | 3 | 2 |  |  |
$4\times 10^{-3}$ | 32 | 5 | 2 | 2 |  |
$2\times 10^{-3}$ | 152 | 26 | 4 | 2 |  |
$10^{-3}$ | 691 | 109 | 18 | 4 | 2 |
$8\times 10^{-4}$ | 841 | 129 | 23 | 3 | 2 |
$6\times 10^{-4}$ | 1610 | 231 | 40 | 8 | 2 |
$4\times 10^{-4}$ | 3791 | 589 | 104 | 15 | 3 |
$2\times 10^{-4}$ | 15859 | 2344 | 375 | 62 | 13 | 2
| Level $\ell$
---|---
$\epsilon$ | 0 | 1 | 2 | 3 | 4 | 5
$8\times 10^{-3}$ | 10 | 2 | 2 |  |  |
$6\times 10^{-3}$ | 19 | 3 | 3 |  |  |
$4\times 10^{-3}$ | 38 | 5 | 4 |  |  |
$2\times 10^{-3}$ | 121 | 18 | 6 | 2 |  |
$10^{-3}$ | 697 | 82 | 49 | 8 |  |
$8\times 10^{-4}$ | 1446 | 118 | 91 | 27 | 6 |
$6\times 10^{-4}$ | 2070 | 218 | 133 | 21 | 3 |
$4\times 10^{-4}$ | 5075 | 484 | 315 | 61 | 14 |
$2\times 10^{-4}$ | 25871 | 1961 | 1668 | 430 | 85 | 14

Table 2: The optimal sample size estimation for MC-FE (left), uniform MLMC-FE
(middle), and adaptive MLMC-FE (right). The simulations were conducted for a
variety of choices of $\epsilon$. The computational cost associated with a
tolerance of $\epsilon=2\times 10^{-4}$ for Monte Carlo was prohibitive; the
entry in the table for this tolerance (with an asterisk) is an estimate.

| MC-FE | Uniform MLMC-FE | Adaptive MLMC-FE
---|---|---|---
$\epsilon$ | Time |  |  Time |  Speedup
---|---
|  Time |  Speedup
---|---
$8\times 10^{-3}$ | 5.67e+00 |  | 4.52e+00 | 1.3
---|---
| 3.50e+00 | 1.6
---|---
$6\times 10^{-3}$ | 7.69e+00 |  | 5.25e+00 | 1.5
---|---
| 5.51e+00 | 1.4
---|---
$4\times 10^{-3}$ | 1.30e+02 |  | 2.32e+01 | 5.6
---|---
| 9.01e+00 | 14.4
---|---
$2\times 10^{-3}$ | 4.83e+02 |  | 4.62e+01 | 10.5
---|---
| 3.16e+01 | 15.3
---|---
$10^{-3}$ | 9.22e+03 |  | 2.47e+02 | 37.3
---|---
| 1.84e+02 | 50.1
---|---
$8\times 10^{-4}$ | 1.50e+04 |  | 2.48e+02 | 60.5
---|---
| 5.73e+02 | 26.2
---|---
$6\times 10^{-4}$ | 2.48e+04 |  | 4.13e+02 | 60.0
---|---
| 6.30e+02 | 39.4
---|---
$4\times 10^{-4}$ | 5.68e+04 |  | 9.29e+02 | 61.1
---|---
| 1.79e+03 | 31.7
---|---
$2\times 10^{-4}$ | 9.62e+05∗ |  | 4.21e+03 | 228.5
---|---
| 1.21e+04 | 79.6
---|---

Table 3: The CPU time in seconds for MC-FE (left), uniform MLMC-FE (middle),
and adaptive MLMC-FE (right), together with speedups for the multilevel
methods, for a variety of choices of $\epsilon$. The computational cost
associated with a tolerance of $\epsilon=2\times 10^{-4}$ for Monte Carlo was
prohibitive; the entry in the table for this tolerance (with an asterisk) is
an estimate.

###  4.3 Properties of geometric parameters

Next, we will explore the plasma boundaries and geometric descriptors of the
expected poloidal flux $\psi$ resulting from the three methods. To ensure a
fair comparison, we will use the results obtained from the MC-FE on the finest
uniform mesh as a reference benchmark.

Plasma boundary. To ascertain the expected location of the plasma boundary, we
first determine the expected solution to the free boundary problem,
$\mathbb{E}[\psi]$, and determine the boundary of this deterministic function.
This boundary is depicted, in dark violet, in Figure 4 along with the plasma
boundaries obtained from 50 random currents, which are shown in light violet
curves. In Figure 5 we present plots depicting the x-points and plasma
boundaries of the expected solution $\psi$ computed using only samples and
corrections from increasingly finer grids for both uniform and adaptive MLMC-
FE approaches. The data was obtained with tolerance $\epsilon=4\times
10^{-4}$. As can be seen when moving from left to right in Figure 5, the
result obtained using the information from the coarsest level (leftmost
column) is progressively corrected with information from increasingly finer
grids, leading to the desired result depicted in the rightmost column.

Among the three methods, MC-FE yields the smoothest plasma boundaries in the
vicinity of the x-point, followed by adaptive MLMC-FE, while the MLMC-FE
approach on geometry-conforming uniform meshes manifests the most pronounced
irregularity in the plasma boundary. The boundary of the expected solution
generated with the uniform grid MLMC-FE method exhibits irregularities as can
be seen in Figure 4. These large deformations can be primarily attributed to
the additional challenges arising from the curved boundaries. We will address
this point in more detail at the end of this section. The top row of Figure 5
demonstrates that using a geometry-conforming mesh provides a more accurate
approximation of the curved structure (in black) of the configuration than
that in the bottom row. These observations underscore the challenge of
striking a balance between preserving geometric fidelity when dealing with
curved boundaries and the desired statistical accuracy of the solution.

|  |  |  |
---|---|---|---|---
|  |  |  |
MC simulation |  |  Uniform MLMC |  |  Adaptive MLMC
Geometry-conforming mesh |  |  Geometry-conforming mesh |  |
Figure 4: The overlayed plasma boundaries of 50 random realizations are displayed in the top row as violet curves. The solid violet line is the plasma boundary of the expected poloidal flux generated with tolerance $\epsilon=4\times 10^{-4}$. The inner and outer walls of the reactor are displayed in solid black and dark red respectively. The bottom row shows the regions close to the x-points in more detail. The dark green dots are the x-points of the expected solution. Each column from left to right corresponds to: simulation with the Monte Carlo approach, MLMC simulation on geometry-conforming uniform meshes, and adaptive MLMC simulation. All simulations were performed using the discretization level $\ell=5$.   Uniform MLMC  |  |  |  |  |
---|---|---|---|---|---
  Adaptive MLMC  |  |  |  |  |
|  $\ell=0$ |  $\ell=1$ |  $\ell=2$ |  $\ell=3$ |  $\ell=4$
Figure 5: The violet curves represent the expected plasma boundaries for
simulations using increasingly finer grids when $\epsilon=4\times 10^{-4}$
using the sample sizes specified in Table 2. Each sub-plot focuses on a region
near the x-point, maintaining the same zoom-in ratio as the second row of
Figure 4. The dark green dots denote the locations of the x-point. The top row
shows the results of MLMC-FE on a set of geometry-conforming uniform meshes,
while the bottom row displays the results for adaptive MLMC.

Geometric descriptors. Table 4 reports some geometric parameters of the
expected poloidal flux $\mathbb{E}[\psi]$ in (3). It is observed that these
parameters are consistent across different simulation techniques, with
agreement typically up to two or, in some cases, one significant digit. Having
been derived from $\mathbb{E}[\psi]$ all these values are deterministic. There
is, however, uncertainty associated with the corresponding quantities derived
from $\psi$, as they are themselves random variables. Uncertainty on these
quantities can be assessed by computing Monte Carlo estimates of their
expectations and variances. The multi-level methodology can be used for this
purpose without any modification by simply regarding them as quantities of
interest in their own right. All the relevant descriptors are computed for
each of the samples gathered during the computational loop at the various
discretization levels and their expectations and variances are estimated
following the same process described on Algorithm 2. The expectations and
variances obtained in this fashion are shown in Table 5.

Despite the advantages of low computational cost, the MLMC-FE-based methods
may encounter difficulties in accurately determining the locations of x-points
and magnetic axis. Note that the x-points, which correspond to saddle points
of the piecewise linear approximation of $\psi$, can only be located at the
nodes of the mesh. The numerical identification of their exact locations,
which often relies on changes in the sign of the discrete gradient, can be
challenging; see [5, 12, 28] for discussions of the computational
difficulties.

In summary, simulations using the uniform MLMC-FE on non-nested geometry-
conforming uniform meshes may encounter a substantial challenge in accurately
identifying the x-point and achieving less accurate quantities, especially for
the plasma boundary, when compared to the results obtained from MC-FE. On the
other hand, the adaptive MLMC-FE approach on a nested adaptively refined mesh
set produces results that closely align with the MC-FE at a much lower
computational cost.

| MC-FE | Uniform MLMC-FE | Adaptive MLMC-FE
---|---|---|---
x point | (5.14,-3.29) | (5.14,-3.29) | (5.14,-3.28)
magnetic axis | (6.41,0.61) | (6.44,0.56) | (6.46,0.54)
strike | (4.16,-3.71) | (4.16,-3.71) | (4.16,-3.71)
points | (5.56,-4.22) | (5.56,-4.22) | (5.56,-4.21)
inverse aspect ratio | 0.32 | 0.32 | 0.32
elongation | 1.86 | 1.87 | 1.86
upper triangularity | 0.43 | 0.43 | 0.43
lower triangularity | 0.53 | 0.53 | 0.53

Table 4: Geometric parameters of the expected poloidal flux $\mathbb{E}[\psi]$
from MC-FE, MLMC-FE with geometry-conforming uniform mesh set, and adaptive
MLMC-FE. The results are generated with an nMSE $4\times 10^{-4}$.

| MC | MLMC
---|---|---
| - | Uniform | Adaptive
Geometric parameters | Mean | Variance | Mean | Variance | Mean | Variance
x point | (5.14,-3.29) | (2.37e-04,1.44e-03) | (5.14,-3.28) | (2.44e-04,6.33e-04) | (5.14,-3.29) | (4.03e-04,2.91e-03)
magnetic axis | (6.41,0.61) | (1.01e-03,6.11e-04) | (6.41,0.60) | (2.65e-02,6.28e-03) | (6.42,0.60) | (3.15e-02,5.81e-03)
strike | (4.17,-3.71) | (4.46e-03,2.23e-03) | (4.16,-3.71) | (1.23e-04,4.30e-05) | (4.17,-3.71) | (4.18e-04,2.00e-04)
points | (5.56,-4.22) | (2.88e-08,3.28e-03) | (5.56,-4.22) | (4.05e-09,5.54e-05) | (5.56,-4.21) | (2.02e-08,2.98e-04)
inverse aspect ratio | 0.32 | 4.76e-06 | 0.32 | 4.23e-06 | 0.32 | 4.31e-06
elongation | 1.86 | 1.50e-04 | 1.86 | 1.83e-03 | 1.86 | 2.29e-03
upper triangularity | 0.43 | 2.49e-04 | 0.43 | 1.45e-03 | 0.43 | 2.08e-03
lower triangularity | 0.53 | 1.08e-04 | 0.53 | 1.52e-03 | 0.53 | 1.31e-03

Table 5: Sample mean of the geometric parameters extracted from each
realization during three simulations MC-FE, uniform MLMC-FE, and adaptive
MLMC-FE. The results are generated with an nMSE $4\times 10^{-4}$.

Meshing curved domains and their effect on MLMC estimations. In the
deterministic setting this geometric error has the undesired consequence of
hindering the decay of the discretization error since, as the mesh is refined,
the discretized computational domain does not converge to the semicircle
bounded by $\Gamma$. In the stochastic setting the geometric error manifests
itself in rendering the Monte Carlo estimator biased and inconsistent. The
inconsistency stems from the fact that, as both the sample size and the mesh
level increase, the Monte Carlo estimator does not converge to the expectation
of the random variable $u$ satisfying the free boundary problem. Instead, the
estimator converges to the expectation of the random variable that satisfies a
perturbation of (28) where the curve $\Gamma$ is not a semicircle, but the
initial polygonal approximation. If the initial mesh is fine enough, this
geometric bias will likely be too small to affect the estimation.

On the other hand, if an exact descriptor of the curved boundaries is
available, the aforementioned difficulty can be overcome by re-sampling the
curved boundaries when building the sequence of finer grids, thus allowing for
a resolution of the curved structures consistent with the respective mesh
parameter. If an exact descriptor is not available it is possible to
approximate it with, for instance, a cubic spline representation that
interpolates the original polygonal representation. This spline surrogate is
then used to re-sample the boundary as the mesh is refined. This strategy,
which gives rise to what we refer to as geometry conforming meshes, was
implemented for the numerical experiments with uniformly refined meshes and
can be seen in use in the top row of Figure 5, where the curved boundary is
represented more accurately as the mesh is refined.

Nevertheless, even if the approximation to the problem geometry is now
consistent with the discretization error, this approach creates additional
challenges. Since the approximations to the curved boundaries are not fixed
across mesh levels, the sequence of meshes is no longer nested—not even in the
case of uniform refinements. Moreover, due to the fact that the sequence of
discrete domains no longer coincides across levels, the domains of definition
of the respective discrete solutions will not overlap, and an extrapolation
step may be needed to compute the multilevel Monte Carlo estimator on a common
computational domain. This strategy, used in our numerical experiments,
introduces an additional extrapolation error. In our case this is evident, for
instance, in the fact that the plasma boundary of the expected solution
$\mathbb{E}[\psi_{h}]$ is considerably less regular in the geometry-conforming
case than it is in the non-geometry-conforming one. This can be seen in Figure
5. The extrapolation error can be taken care of through careful post-
processing. One option is to project or interpolate the numerical solutions
into a subdomain common to all grids so that no extrapolation is needed for
evaluation. This strategy was employed to produce Figure 6 successfully
eliminating the spurious oscillations in the plasma boundary. However, doing
this requires considerable computational work and reduces the time savings
obtained from MLMC.

One further difficulty is that the re-sampling of the boundaries is impossible
to perform in a straightforward fashion in the case of adaptively refined
meshes. Thus, the geometric approximation remains fixed at the initial level
of refinement. This can be seen in the bottom row of Figure 5, where the solid
black line represents the polygonal approximation to the curved boundary of
the divertor. The approximation improves as the mesh is refined for the
uniformly refined mesh, but stays fixed for the adaptive strategy.

|  |  |  |
---|---|---|---|---
$\ell=0$ |  $\ell=1$ |  $\ell=2$ |  $\ell=3$ |  $\ell=4$
Figure 6: The violet curves represent the expected plasma boundaries for post-
processed simulations using increasingly finer grids when $\epsilon=4\times
10^{-4}$ and sample size specified in Table 2. Each sub-plot focuses on a
region near the x-point, maintaining the same zoom-in ratio as the second row
of Figure 4. The dark green dots denote the locations of the x-point.

##  5 Concluding remarks

The objective of this study is to evaluate the performance of MC-FE and
several variants of MLMC-FE for the Grad-Shafranov free boundary problem under
high-dimensional uncertainties in the currents. The plasma equilibrium problem
presents challenges and complexities associated with the physical properties
of the system, the values of the parameters of the experiment, and even the
approximation of geometric features of the model. Examples include the
degradation of the stability of the equilibrium associated with an increase of
plasma elongation or beta; the presence of finite current near the plasma
boundary or, as in this study, the complex boundary structures, which must be
treated carefully in the context of multilevel methods. These physically
grounded challenges typically get reflected in the computations, for instance,
as an increase in the difficulty of a linearized algorithm to converge to a
desired tolerance or as the need for higher precision in the numerical
computation required to resolve the desired features. In other words,
complicated physical scenarios typically result in an increase in
computational costs which might easily become overwhelming when combined with
uncertainty quantification efforts. In cases like these, an MLMC strategy can
help to mitigate the costs required to gather statistical information at the
expense of potentially overlooking subtle behavior that may not be resolved in
the samples gathered in coarse levels of refinement.

It is important to remark that the multilevel Monte Carlo techniques presented
rely only on the independence between different realizations across
discretization levels and not on the particular probability distribution
obeyed by the stochastic parameters. In this study, the current intensities
were modeled as independent uniformly distributed random variables, but the
method would have worked without any modification even if the current
intensities had obeyed correlated distributions as long as realizations across
discretization levels are independent of each other. Moreover, further
stochastic components—such as the parameters appearing in the profile models
in equation (1b)—can be included and handled simultaneously.

In this study, we handled the geometric difficulties in the uniform MLMC-FE
case using a sequence of uniform meshes that conform to the curved boundaries
in the geometry. We also found that a feature of adaptive MLMC-FE is that
adaptive gridding leads to a somewhat better representation of geometric
quantities; in contrast, errors introduced by computations on uniform coarse
meshes lead to some distortions of some features. The traditional advantage of
adaptive refinement reduced computational cost, was less clearly present. But
the overall advantage of the MLMC approach in reducing costs is dramatic.

##  6 Acknowledgements

All the free boundary computations were done using the code FEEQS.M [23, 24].
The authors are deeply grateful to Holger Heumann, the INRIA CASTOR team, and
all the development team of the CEDRES++ free boundary solver for sharing
access to the code and for helping us get up to speed with its usage.

Howard Elman has been partially supported by the U. S. Department of Energy
under grant DE-SC0018149 and by the U. S. National Science Foundation under
grant DMS1819115. Jiaxing Liang has been partially funded by the U. S.
Department of Energy under grant DE‐SC0018149. Tonatiuh Sánchez-Vizuet has
been partially funded by the U. S. National Science Foundation through the
grant NSF-DMS-2137305.

##  7 Appendix

We now describe some technical mathematical and implementation details related
to the work presented. This additional material is provided for the sake of
completeness but is not needed to follow the discussion from the previous
sections.

###  7.1 Weak formulation of the Grad-Shafranov equation

Deterministic problem. We will start by introducing the weak formulation for
the Grad-Shafranov problem in the case where all the parameters involved are
deterministic. Consider a semi-circle centered at the origin, boundary
$\Gamma$ and with radius $\rho$ such that it fully contains all the relevant
reactor components depicted in Figure 1. If $\Omega$ denotes the region
surrounded by $\Gamma$ then, by construction, for any in $\Omega^{c}$ the
right-hand side of (1a) will vanish identically. We will then, following [22],
consider the space of real-valued functions

$Z:=\left\\{u:\Omega\rightarrow\mathbb{R}\,\Bigg{|}\,\int_{\Omega}u^{2}x\,dxdy<\infty\,,\,\int_{\Omega}\frac{\left|\nabla
u\right|^{2}}{x}dxdy<\infty;\,\text{ and }\,u(0,y)=0\right\\}\cap
C^{0}\left(\overline{\Omega}\right).$ (26)

This space arises naturally when testing equation (1a) using the weighted
$L_{2}$ inner product defined by

$\langle u,v\rangle:=\int_{\Omega}uvx\,dxdy,$

which leads to the finite energy requirement appearing in the second
inequality in the definition (26). The third requirement ($u(0,y)=0$) is a
result of the anti-symmetry of the problem with respect to reflections across
the axis of symmetry of the reactor and has the effect of ensuring that the
quantity

$\left\|u\right\|_{Z}:=\left(\int_{\Omega}\frac{\left|\nabla
u\right|^{2}}{x}dxdy\right)^{1/2}$ (27)

does indeed define a norm in the space $Z$. We will refer to this norm as the
energy norm. The space $Z$ defined in (26) endowed with the energy norm (27)
is the natural function space to look for variational solutions to the
deterministic linearized Grad-Shafranov equation (1). The variational
formulation of the problem, as derived in [1, 17, 29], is that of finding
$\psi\in Z$ such that for every test $\varphi\in Z$ it holds that:

$\displaystyle\int_{\Omega_{p}}\frac{1}{\mu
x}\nabla\psi\cdot\nabla\varphi\,dxdy-\int_{\Omega_{p}}\left(rp^{\prime}(\psi)+\frac{1}{\mu_{0}x}ff^{\prime}(\psi)\right)\varphi\,dxdy+\int_{\Gamma}\psi\,N\,\varphi\,dc$
$\displaystyle\hskip
14.45377pt+\int_{\Gamma}\int_{\Gamma}\left(\psi(\boldsymbol{x}_{1})-\psi(\boldsymbol{x}_{2})\right)M(\boldsymbol{x}_{1},\boldsymbol{x}_{2})\left(\varphi(\boldsymbol{x}_{1})-\varphi(\boldsymbol{x}_{2})\right)\,dc(\boldsymbol{x}_{1})\,dc(\boldsymbol{x}_{2})=\sum_{k=1}^{M_{c}}\frac{I_{k}}{S_{k}}\int_{\Omega_{C_{k}}}\varphi\,dxdy.$
(28)

Above, the magnetic permeability $\mu$ is either a function of $\psi$ inside a
region occupied by a ferromagnetic material,
$\mu=\mu(|\nabla\psi|^{2}/x^{2})$, or a constant $\mu=\mu_{0}$ elsewhere,
$\boldsymbol{x}_{i}=(x_{i},y_{i})$ denote position vectors, $\Omega_{C_{k}}$
denotes the region occupied by the $k$-th coil, and the total number of coils
is denoted by $M_{c}$. The following quantities, appearing on the left-hand
side of (28), are related to the Green’s function associated with the operator
on the left-hand side of (1a)

$\displaystyle
N:=\frac{1}{x}\left(\frac{1}{\delta_{+}}+\frac{1}{\delta_{-}}-\frac{1}{\rho}\right)$
$\displaystyle\,,\qquad\delta_{\pm}:=\sqrt{x^{2}+(\rho\pm
y)^{2}}\,,\qquad\kappa(\boldsymbol{x}_{1},\boldsymbol{x}_{2}):=\sqrt{\frac{4x_{1}x_{2}}{(x_{1}+x_{2})^{2}+(y_{1}-y_{2})^{2}}}\,,$
$\displaystyle M(\boldsymbol{x}_{1},\boldsymbol{x}_{2}):=\,$
$\displaystyle\frac{\kappa(\boldsymbol{x}_{1},\boldsymbol{x}_{2})}{2\pi(x_{1}x_{2})^{3/2}}\left(\frac{2-\kappa^{2}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})}{2-2\kappa^{2}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})}E(\kappa(\boldsymbol{x}_{1},\boldsymbol{x}_{2}))-K(\kappa(\boldsymbol{x}_{1},\boldsymbol{x}_{2}))\right)\,,$

where $E(\kappa(\boldsymbol{x}_{1},\boldsymbol{x}_{2}))$ and
$K(\kappa(\boldsymbol{x}_{1},\boldsymbol{x}_{2}))$ are complete elliptic
integrals of the first and second kind respectively [30].

Accounting for stochasticity. We now consider the stochasticity in the
currents and allow the vector of currents to be a $d$-dimensional random
variable $\boldsymbol{\omega}$ uniformly distributed over the parameter space
$W$ defined in (2). It is clear that in this case for any particular
realization of the currents $\boldsymbol{\omega}$ we will obtain a different
equilibrium configuration $\psi(\boldsymbol{\omega})$ that belongs to the
Banach space $Z$ defined in (26). Moreover, since for every
$\boldsymbol{\omega}\in W$ the resulting equilibrium configuration has finite
energy, it then holds for the expected value of the equilibrium that

$\mathbb{E}\left(\|\psi\|_{Z}^{2}\right)<\infty.$

In mathematical terms, we say that the stochasticity of the currents
transforms the solution $\psi$ to (1) into a Banach space-valued random
variable with finite expected energy. If $\boldsymbol{\omega}$ belongs to a
complete and separable probability space $\left(W,\Sigma,\mathbb{P}\right)$,
the class of such random variables forms what is known as a Bochner space [9].
In our particular case, solutions to (1) are mappings from the parameter space
$W$ to the Banach space $Z$ that, as functions of $\boldsymbol{\omega}$,
belong to the Bochner space

$L^{2}(W,\Sigma,\mathbb{P};Z):=\\{u:W\rightarrow Z\;\big{|}\;u\text{ strongly
measurable,}\;\;\|u\|_{L^{2}(W,Z)}<\infty\\},$

where the norm $\|\cdot\|_{L^{2}(W,Z)}$ is precisely defined in terms of the
expected energy

$\|u\|_{L^{2}(W,Z)}:=\left(\int_{W}\left\|u(\cdot,\omega)\right\|_{Z}^{2}d\mathbb{P}(\omega)\right)^{1/2}=\left(\mathbb{E}\left(\left\|u(\cdot,\omega)\right\|_{Z}^{2}\right)\right)^{1/2}.$

###  7.2 Mathematical results on multilevel Monte Carlo

The proceeding analysis on the computational cost to compute the MLMC-FE
estimator in terms of the desired relative accuracy $\epsilon$ is based on the
following Theorem (established and proven in [8]), which quantifies the
distribution of the computational effort across discretization levels in terms
of the relation between the decay rate of the variance of the correction
terms, $V_{\ell}$, and the increase of the computational cost, $C_{\ell}$, as
the mesh is refined. In particular, it states that if the variance of the
correction terms decays faster than the increase of computational cost, the
dominant computational expense takes place on the coarsest grid—see also [19].

######  Theorem 1.

Suppose there exist positive constants $a,b,c$ such that
$a\geq\frac{1}{2}\min(b,c)$,

  1. (i)

$\left\|\mathbb{E}\left(u-u_{\ell}\right)\right\|_{Z}=\mathcal{O}\left(M_{\ell}^{-a}\right)$,

  2. (ii)

$V_{\ell}=\mathcal{O}\left(M_{\ell}^{-b}\right)$,

  3. (iii)

$C_{\ell}=\mathcal{O}\left(M_{\ell}^{c}\right)$.

Then for any positive $\epsilon<e^{-1}$ small enough, there exists level $L$
and sample size $N_{\ell}$ for which the multilevel estimator
$A_{\text{MLMC}}$ has an nMSE with

$\frac{\mathbb{E}\left[\left\|\mathbb{E}(u)-A_{\text{MLMC}}(u_{L})\right\|_{Z}^{2}\right]}{\mathbb{E}\left[\left\|\mathbb{E}(u)\right\|_{Z}^{2}\right]}<\epsilon^{2},$

and the total computation cost $C$ with bound

$C\left(A_{\text{MLMC}}\right)=\left\\{\begin{array}[]{ll}\mathcal{O}\left(\epsilon^{-2}\right),&b>c,\\\
\mathcal{O}\left(\epsilon^{-2}\left(\log\epsilon\right)^{2}\right),&b=c,\\\
\mathcal{O}\left(\epsilon^{-2-\frac{(c-b)}{a}}\right),&0<b<c.\end{array}\right.$

## References

  * [1] R. Albanese, J. Blum, and O. D. Barbieri.  On the solution of the magnetic flux equation in an infinite domain.  In S. Methfessel and G. Thomas, editors, 8th. European conference on Computational Physics: Computing in plasma physics, volume 10D of Europhysics Conference Abstracts, pages 41–44. European Physical Society, 1986.
  * [2] V. Amoskov, A. Belov, V. Belyakov, Y. Gribov, A. Kavin, V. Kukhtin, E. Lamzin, K. Lobanov, N. Maximenkova, A. Mineev, and S. Sytchevsky.  Stray magnetic field at plasma initiation produced by ferromagnetic elements of the ITER tokamak complex.  Plasma Devices and Operations, 17(4):238–249, Dec. 2009.
  * [3] I. Babuška, F. Nobile, and R. Tempone.  A stochastic collocation method for elliptic partial differential equations with random input data.  SIAM Journal on Numerical Analysis, 45(3):1005–1034, 2007.
  * [4] A. Barth, C. Schwab, and N. Zollinger.  Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients.  Numerische Mathematik, 119(1):123–161, 2011.
  * [5] C. Bierig and A. Chernov.  Estimation of arbitrary order central statistical moments by the multilevel Monte Carlo method.  Stochastic Partial Differential Equations. Analysis and Computations, 4(1):3–40, 2016.
  * [6] C. Carstensen.  An adaptive mesh-refining algorithm allowing for an h 1 stable l 2 projection onto courant finite element spaces.  Constructive Approximation, 20(4):549–564, Oct. 2003.
  * [7] C. Carstensen, M. Feischl, M. Page, and D. Praetorius.  Axioms of adaptivity.  Computers & Mathematics with Applications, 67(6):1195–1253, Apr. 2014.
  * [8] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup.  Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients.  Computing and Visualization in Science, 14(1):3–15, 2011.
  * [9] G. Da Prato and J. Zabczyk.  Encyclopedia of Mathematics and its Applications: Stochastic Equations in Infinite Dimensions, Series Number 152.  Cambridge University Press, Cambridge, England, 2 edition, Apr. 2014.
  * [10] W. Dörfler.  A convergent adaptive algorithm for Poisson’s equation.  SIAM Journal on Numerical Analysis, 33(3):1106–1124, 1996.
  * [11] M. Eigel, C. Merdon, and J. Neumann.  An adaptive multilevel Monte Carlo method with stochastic bounds for quantities of interest with uncertain data.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):1219–1245, 2016.
  * [12] H. C. Elman, J. Liang, and T. Sánchez-Vizuet.  Surrogate approximation of the Grad-Shafranov free boundary problem via stochastic collocation on sparse grids.  Journal of Computational Physics, 448:Paper No. 110699, 20, 2022.
  * [13] H. R. Fairbanks, A. Doostan, C. Ketelsen, and G. Iaccarino.  A low-rank control variate for multilevel Monte Carlo simulation of high-dimensional uncertain systems.  Journal of Computational Physics, 341:121–139, 2017.
  * [14] B. Faugeras and H. Heumann.  FEM-BEM coupling methods for Tokamak plasma axisymmetric free-boundary equilibrium computations in unbounded domains.  Journal of Computational Physics, 343:201 – 216, 2017.
  * [15] S. Funken, D. Praetorius, and P. Wissgott.  Efficient implementation of adaptive P1-FEM in Matlab (extended preprint).  Technical report, Institute for Analysis and Scientific Computing, Vienna University of Technology, 2008.
  * [16] S. Funken, D. Praetorius, and P. Wissgott.  Efficient implementation of adaptive P1-FEM in Matlab.  Computational Methods in Applied Mathematics, 11(4):460–490, 2011.
  * [17] G. Gatica and G. Hsiao.  The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems.  Journal of Mathematical Analysis and Applications, 189(2):442 – 461, 1995.
  * [18] R. Ghanem and P. Spanos.  A stochastic Galerkin expansion for nonlinear random vibration analysis.  Probabilistic Engineering Mechanics, 8(3-4):255–264, 1993.
  * [19] M. B. Giles.  Multilevel Monte Carlo methods.  Acta Numerica, 24:259–328, 2015.
  * [20] P.-A. Gourdain, J.-N. Leboeuf, and R. Neches.  High-resolution magnetohydrodynamic equilibrium code for unity beta plasmas.  J. Comput. Phys., 216(1):275–299, July 2006.
  * [21] H. Grad and H. Rubin.  Hydromagnetic equilibria and force-free fields.  In Proc. Second international conference on the peaceful uses of atomic energy, Geneva, volume 31,190, New York, Oct 1958. United Nations.
  * [22] V. Grandgirard.  Modelisation de l’equilibre d’un plasma de tokamak.  Technical report, Université de Franche-Comté, 1999.
  * [23] H. Heumann.  FEEQS.M.  http://www-sop.inria.fr/members/Holger.Heumann/Software.html.
  * [24] H. Heumann, J. Blum, C. Boulbe, B. Faugeras, G. Selig, J.-M. Ané, S. Brémond, V. Grandgirard, P. Hertout, E. Nardon, and et al.  Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications.  Journal of Plasma Physics, 81(3):905810301, 2015.
  * [25] H. Heumann and F. Rapetti.  A finite element method with overlapping meshes for free-boundary axisymmetric plasma equilibria in realistic geometries.  Journal of Computational Physics, 334:522 – 540, 2017.
  * [26] H. Hoel, E. von Schwerin, A. Szepessy, and R. Tempone.  Adaptive multilevel Monte Carlo simulation.  In Numerical Analysis of Multiscale Computations, volume 82 of Lect. Notes Comput. Sci. Eng., pages 217–234. Springer, Heidelberg, 2012.
  * [27] H. Hoel, E. von Schwerin, A. Szepessy, and R. Tempone.  Implementation and analysis of an adaptive multilevel Monte Carlo algorithm.  Monte Carlo Methods and Applications, 20(1):1–41, 2014.
  * [28] F. Hofmann.  FBT - a free-boundary tokamak equilibrium code for highly elongated and shaped plasmas.  Comput. Phys. Commun., 48(2):207–221, Feb. 1988.
  * [29] G. C. Hsiao and S. Zhang.  Optimal order multigrid methods for solving exterior boundary value problems.  SIAM Journal on Numerical Analysis, 31(3):680–694, 1994.
  * [30] J. D. Jackson.  Classical Electrodynamics.  Wiley, New York, NY, 2nd edition, 1975.
  * [31] J. Johnson, H. Dalhed, J. Greene, R. Grimm, Y. Hsieh, S. Jardin, J. Manickam, M. Okabayashi, R. Storer, A. Todd, D. Voss, and K. Weimer.  Numerical determination of axisymmetric toroidal magnetohydrodynamic equilibria.  Journal of Computational Physics, 32(2):212–234, Aug. 1979.
  * [32] A. Khodadadian, M. Parvizi, and C. Heitzinger.  An adaptive multilevel Monte Carlo algorithm for the stochastic drift-diffusion-Poisson system.  Computer Methods in Applied Mechanics and Engineering, 368:113163, 23, 2020.
  * [33] R. Kornhuber and E. Youett.  Adaptive multilevel Monte Carlo methods for stochastic variational inequalities.  SIAM Journal on Numerical Analysis, 56(4):1987–2007, 2018.
  * [34] S. Larsson and V. Thomée.  Partial Differential Equations with Numerical Methods, volume 45 of Texts in Applied Mathematics.  Springer-Verlag, Berlin, 2003.
  * [35] F. Law, A. Cerfon, and B. Pehersdorfer.  Accelerating the estimation of collisionless energetic particle confinement statistics in stellarators using multifidelity Monte Carlo.  Nucl. Fusion, 62:076019, 2022.
  * [36] R. Lüst and A. Schlüter.  Axialsymmetrische magnetohydrodynamische Gleichgewichtskonfigurationen.  Z. Naturf, 12a:850–854, 1957.
  * [37] J. Luxon and B. Brown.  Magnetic analysis of non-circular cross-section tokamaks.  Nuclear Fusion, 22(6):813–821, jun 1982.
  * [38] N. Metropolis and S. Ulam.  The Monte Carlo method.  Journal of the American Statistical Association, 44:335–341, 1949.
  * [39] P. Morin, R. H. Nochetto, and K. G. Siebert.  Data oscillation and convergence of adaptive FEM.  SIAM Journal on Numerical Analysis, 38(2):466–488, Jan. 2000.
  * [40] M. Motamed and D. Appelö.  A multiorder discontinuous Galerkin Monte Carlo method for hyperbolic problems with stochastic parameters.  SIAM Journal on Numerical Analysis, 56(1):448–468, 2018.
  * [41] F. Nobile and F. Tesei.  A multi level Monte Carlo method with control variate for elliptic PDEs with log-normal coefficients.  Stochastic Partial Differential Equations. Analysis and Computations, 3(3):398–444, 2015.
  * [42] V. D. Shafranov.  On magnetohydrodynamical equilibrium configurations.  Soviet Physics JETP, 6:545–554, 1958.
  * [43] J. R. Shewchuk.  Delaunay refinement algorithms for triangular mesh generation.  Computational Geometry, 22(1-3):21–74, May 2002.
  * [44] A. L. Teckentrup, R. Scheichl, M. B. Giles, and E. Ullmann.  Further analysis of multilevel monte carlo methods for elliptic PDEs with random coefficients.  Numerische Mathematik, 125(3):569–600, Mar. 2013.
  * [45] B. P. Welford.  Note on a method for calculating corrected sums of squares and products.  Technometrics, 4(3):419–420, 1962.
  * [46] D. Xiu.  Numerical Methods for Stochastic Computations.  Princeton University Press, Princeton, NJ, 2010.  A spectral method approach.