2. Scattering geometry and ScaRF code

The coordinate system ScaRF code is based on the toroidal plasma configuration of figure 1, where two separate coordinate systems (CSs) are present relative to the toroidal plasma configuration. The $(x_B, y_B, z_B)$ CS corresponds to the magnetic flux density CS where the $z_B$ direction corresponds to the direction of the magnetic field. The $(x_p, y_p, z_p)$ CS corresponds to the plasma coordinate system. The magnetic flux density can have a poloidal and toroidal component therefore the two CSs can be related via Euler-rotation angles, as shown in figure 2. Following the work of Stix (Reference Stix1992)), the relative permittivity of the cold plasma is given by the following equation in the $(x_B, y_B, z_B)$ CS:

(2.1)\begin{equation} {{\tilde{\boldsymbol \varepsilon}}}_B = \left[\begin{array}{@{}ccc@{}} S & -\textrm{i} D & 0 \\ +\textrm{i} D & S & 0 \\ 0 & 0 & P \end{array}\right], \end{equation}

where the $S$, $D$, $P$ parameters are defined as

(2.2)\begin{equation} \left.\begin{gathered} S= \tfrac{1}{2} \left( R + L \right),\quad D = \tfrac{1}{2} \left( R - L \right), \\ P= 1 + P_e + P_i, \quad R = 1 + R_e + R_i, \\ L= 1 + L_e + L_i, \quad R_e =\frac{1 +i n_{\textrm{col}}}{1 + C_e + i n_{\textrm{col}}} P_e,\\ R_i =\frac{1 +i n_{\textrm{col}}}{1 + C_i + i n_{\textrm{col}}} P_i, \quad L_e = \frac{1 +i n_{\textrm{col}}}{1 - C_e + i n_{\textrm{col}}} P_e, \\ L_i= \frac{1 +i n_{\textrm{col}}}{1 - C_i + i n_{\textrm{col}}} P_i,\quad P_e ={-}\left(\frac{\omega_{\textrm{pe}}}{\omega}\right)^{2} \frac{1}{1+i n_{\textrm{col}}}, \\ P_i={-}\left(\frac{\omega_{\textrm{pi}}}{\omega}\right)^{2} \frac{1}{1+i n_{\textrm{col}}}, \quad C_e ={-}\frac{\omega_{\textrm{ce}}}{\omega}, \\ C_i= \frac{\omega_{\textrm{ci}}}{\omega}, \quad \omega_{\textrm{pe}} = q_e \sqrt{ \frac{n_e}{m_e \varepsilon_0} }, \\ \omega_{\textrm{pi}}= q_i \sqrt{ \frac{n_i}{m_i \varepsilon_0} } = Zq_e \sqrt{ \frac{n_i}{A m_p \varepsilon_0} }, \quad \omega_{\textrm{ce}} = \frac{q_e B_{z_B} }{m_e}, \\ \omega_{\textrm{ci}}= \frac{ Z q_e B_{z_B} }{A m_p}, \quad n_{\textrm{col}} = \frac{\nu_{\textrm{col}}}{\omega}. \end{gathered}\right\} \end{equation}

The variables in the above equations are: the magnitude of the electron charge $q_e$, the electron rest mass $m_e$, the atomic number $Z$, the atomic mass number $A$, the proton rest mass $m_p$, the permittivity of free space $\varepsilon _0$, the frequency of the electromagnetic radiation $\omega = 2{\rm \pi} f$ (where $f$ is the frequency in Hz), the electron and ion plasma densities $n_e$ and $n_i$, respectively (in $\textrm {m}^{-3}$), and the electron and ion collision rate $\nu _{\textrm {col}}$. The frequencies $\omega _{\textrm {pe}}$ and $\omega _{\textrm {pi}}$ are the electron and ion plasma resonant frequencies and $\omega _{\textrm {ce}}$ and $\omega _{\textrm {ci}}$ are the electron and ion cyclotron frequencies, respectively. In this work $\omega _{\textrm {ce}}=7.915\times 10^{11}\ \textrm {rad}\ \textrm {s}^{-1}$ and $\omega _{\textrm {ci}}=2.155\times 10^{8}\ \textrm {rad}\ \textrm {s}^{-1}$. Collisional absorption is absent in the current model, but can be easily included by changing the relative permittivity tensor in (2.1) accordingly. The relative permittivity tensor of (2.1) can be expressed in the plasma CS $(x_p, y_p, z_p)$ by use of the transformation:

(2.3)\begin{equation} \boldsymbol{\tilde{\boldsymbol \varepsilon}}_p = {\boldsymbol{\tilde{{\boldsymbol{\mathsf{M}}}}}}^{{-}1} {\boldsymbol{\tilde{\boldsymbol \varepsilon}}}_B {\boldsymbol{\tilde{{\boldsymbol{\mathsf{M}}}}}}, \end{equation}

where ${\boldsymbol{\tilde {\boldsymbol \varepsilon }}}_B$ is defined in (2.1) and $\boldsymbol{\tilde {{\boldsymbol{\mathsf{M}}}}}$ is the Euler-angle transformation matrix. If a ripple is present at the torus–plasma surface, diffraction of the incident microwave radiation could occur. To study this effect, a periodic ripple at the plasma torus surface is considered. Its geometry is shown in figure 3, where the ripple has been defined as

(2.4)\begin{equation} h(x) = d - d \cos\left( \frac{2{\rm \pi}}{\varLambda} x \right), \end{equation}

where $h(x)$ is the cosinusoidally varying ripple height, $d$ is the ripple amplitude and $\varLambda$ the spatial period of the ripple (spatial mode).

Figure 1. A tokamak plasma torus is shown with the corresponding coordinate systems of interest. The coordinate system $(x_B, y_B, z_B)$ corresponds to the magnetic flux density, $\boldsymbol {B}$, coordinate system while $(x_p, y_p, z_p)$ corresponds to the plasma coordinate system. The $z_p$-component of the magnetic flux density corresponds to the toroidal magnetic flux density component, $B_{\textrm {tor}}$, while the $x_py_p$-component corresponds to the poloidal magnetic flux density component, $B_{\textrm {pol}}$ (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019).

Figure 2. The relation between the $(x_B, y_B, z_B)$ and the $(x_p, y_p, z_p)$ coordinate systems. The angles $\varphi _B$, $\theta _B$ and $\psi _B$ are the Euler angles that connect the two coordinate systems. All the angles are defined positive as counter-clockwise (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019).

Figure 3. A plasma ripple at the torus boundary is considered as a periodic spatial modulation, i.e. as a plasma grating. The microwave radiation is represented as a plane wave incident from the top towards the bottom region. The incident wavevector is shown as ${\boldsymbol {k}}_{\textrm {inc}}$ and the incident angles are defined as $\phi$ and $\theta$. The scattering CS $(x,y,z)$ is related to the plasma CS by $x=y_p$, $y=z_p$ and $z=x_p$, as shown in the figure. The plasma relative permittivities of the top and of the bottom regions are defined as ${\boldsymbol{\tilde {\boldsymbol \varepsilon }}}_1$ and ${\boldsymbol{\tilde {\boldsymbol \varepsilon }}}_2$, respectively. The plasma grating is assumed to have a sinusoidal profile of periodicity $\varLambda$ along the $x$ direction, and an amplitude spatial variation of $d$ (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019).

The incident microwave excitation is considered to be a plane wave with an azimuthal angle of incidence $\phi$ and a polar angle of incidence $\theta$ (figure 3). The regions above and below the periodic interface are plasma regions of different plasma densities. Their tensor relative permittivities ${\boldsymbol{\tilde{\boldsymbol \varepsilon}}}_{l,p}$ can be determined by suitable use of (2.3) in the plasma CS. Then the relative permittivities ${\boldsymbol{\tilde {\boldsymbol \varepsilon }}}_1$ and ${\boldsymbol{\tilde {\boldsymbol \varepsilon }}}_2$ needed for the diffraction analysis in the $(x,y,z)$ CS can be found from

(2.5)\begin{equation} {\boldsymbol{\tilde{\boldsymbol \varepsilon}}}_{\ell} = {\boldsymbol{ \tilde{ {\boldsymbol{\mathsf{Q }}} } } } {\boldsymbol{\tilde{\boldsymbol \varepsilon}}}_{\ell,p} {\boldsymbol{ \tilde{\boldsymbol{\mathsf{Q}}} } }^\textrm{T}, \quad \ell=1,2, \end{equation}

where $\boldsymbol{\tilde{ {\boldsymbol{\mathsf{Q}}}}}$ is the transformation matrix from the plasma CS $(x_p, y_p, z_p)$ to the $(x,y,z)$ CS of the ScaRF code. The FDFD method in ScaRF is used in the analysis of the plasma grating in figure 3, with FBPBC in the $xy$ plane and PML in the $xy$ planes parallel to the $z$ axis. The FDFD method solves Maxwell's equations in the frequency domain. It is a stable and rigorous method of well-known error sources (Smith Reference Smith1996; Sadiku Reference Sadiku2001; Taflove & Hagness Reference Taflove and Hagness2005), it can model systems of complex geometry and can run in parallel for computationally intensive problems. In the FDFD method, electric and magnetic fields are staggered in space according to the Yee grid (figure 4) and finite differences are used to approximate Maxwell's equations. This results in a large linear algebraic system whose solution provides the electromagnetic fields in space. In particular, after normalization of the magnetic field according to $\boldsymbol {\tilde {H}} \equiv -\textrm {i}Z_{0}\boldsymbol {H}$, where $Z_{0}$ is the free space impedance and $i$ is $\sqrt {-1}$, Maxwell's equations with the wave-absorbing PML layer truncating the computational grid become:

(2.6)\begin{gather} \boldsymbol{\nabla} \times \boldsymbol{E}=k_{0}\left[ \boldsymbol{\tilde{\boldsymbol \varepsilon}} \right]\boldsymbol{\tilde{H}} \end{gather}

(2.7)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{\tilde{H}}=k_{0}\left[\boldsymbol{ \tilde{\mu} }\right]\boldsymbol{E}, \end{gather}

where $[\boldsymbol{ \tilde {\varepsilon }} ] \equiv {\boldsymbol{\mathsf{J}}} [\boldsymbol{ \varepsilon} ] {\boldsymbol{\mathsf{J}}}^{T}/{\mbox {det} ({\boldsymbol{\mathsf{J}}})}$ and $[\boldsymbol{ \tilde {\mu }} ] \equiv {\boldsymbol{\mathsf{J}}}[ \boldsymbol{\mu} ] {\boldsymbol{\mathsf{J}}}^{T}/{\mbox {det} ({\boldsymbol{\mathsf{J}}})}$ are relative permittivity and permeability tensors, defined so as to implement the PML for anisotropic media (Oskooi & Johnson Reference Oskooi and Johnson2011), where $\boldsymbol{\boldsymbol{\mathsf{J}}}\equiv \mbox {diag}(s_{x}^{-1},s_{y}^{-1},s_{z}^{-1})$ and $s_{w}\equiv \kappa _{w}+\textrm {i}({\sigma _{w}}/{\omega })$, $w=\{x,y,z\}$ are the PML stretching factors with $\kappa _{w}\geqslant 1$ as the evanescent wave absorption parameter and $\sigma _{w}$ as the PML conductivity. The parameters of the stretching factors are polynomials (Oskooi & Johnson Reference Oskooi and Johnson2011) spatially varying along the $w$ direction. It is numerically convenient to simplify Maxwell's equations (2.6)–(2.7) by normalization of the grid coordinates as $\tilde {w}=k_{0}w$, $w=\{x,y,z\}$, which leads to

(2.8)\begin{gather} \frac{\partial E_{z}}{\partial \tilde{y}}-\frac{\partial E_{y}}{\partial \tilde{z}}=\tilde{\mu}_{xx}\tilde{H}_{x}+\tilde{\mu}_{xy}\tilde{H}_{y}+\tilde{\mu}_{xz}\tilde{H}_{z} \end{gather}

(2.9)\begin{gather}\frac{\partial E_{x}}{\partial \tilde{z}}-\frac{\partial E_{z}}{\partial \tilde{x}}=\tilde{\mu}_{yx}\tilde{H}_{x}+\tilde{\mu}_{yy}\tilde{H}_{y}+\tilde{\mu}_{yz}\tilde{H}_{z} \end{gather}

(2.10)\begin{gather}\frac{\partial E_{y}}{\partial \tilde{x}}-\frac{\partial E_{x}}{\partial \tilde{y}}=\tilde{\mu}_{zx}\tilde{H}_{x}+\tilde{\mu}_{zy}\tilde{H}_{y}+\tilde{\mu}_{zz}\tilde{H}_{z} \end{gather}

(2.11)\begin{gather}\frac{\partial \tilde{H}_{z}}{\partial \tilde{y}}-\frac{\partial \tilde{H}_{y}}{\partial \tilde{z}}=\tilde{\varepsilon}_{xx}E_{x}+\tilde{\varepsilon}_{xy}E_{y}+\tilde{\varepsilon}_{xz}E_{z} \end{gather}

(2.12)\begin{gather}\frac{\partial \tilde{H}_{x}}{\partial \tilde{z}}-\frac{\partial \tilde{H}_{z}}{\partial \tilde{x}}=\tilde{\varepsilon}_{yx}E_{x}+\tilde{\varepsilon}_{yy}E_{y}+\tilde{\varepsilon}_{yz}E_{z} \end{gather}

(2.13)\begin{gather}\frac{\partial \tilde{H}_{y}}{\partial \tilde{x}}-\frac{\partial \tilde{H}_{x}}{\partial \tilde{y}}=\tilde{\varepsilon}_{zx}E_{x}+\tilde{\varepsilon}_{zy}E_{y}+\tilde{\varepsilon}_{zz}E_{z}. \end{gather}

In ScaRF, the discretization of (2.8)–(2.13) is applied by the approximation of the spatial derivatives by staggered central differences (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019). Then the electric and magnetic fields can be obtained everywhere in the computational domain by solving the linear system $\boldsymbol{\tilde{\boldsymbol{\mathsf{A}}}}[ \boldsymbol {e}, \boldsymbol {\tilde {h}} ]^{\textrm {T}}=\boldsymbol {b}$, where $\boldsymbol {e} \equiv [ \boldsymbol{e_{x}},\boldsymbol{e_{y}},\boldsymbol{e_{x}}]^{\textrm {T}}$, $\boldsymbol {\tilde {h}} \equiv [ \boldsymbol{\tilde {h}_{x}},\boldsymbol{\tilde {h}_{y}},\boldsymbol{\tilde {h}_{x}}]^{\textrm {T}}$ and the right-hand side $\boldsymbol {b}$ is specified from the plane-wave incident source by use of the total field/scattered field (TFSF) technique (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019). The incident source field must satisfy the dispersion and polarization of waves in the cold plasma medium and is analytically described by Papadopoulos et al. (Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019). It is noted that similarly, the right-hand side $\boldsymbol {b}$ can by defined from a superposition of plane waves with amplitudes specified from spectral analysis to represent a spatially localized incident beam.

Figure 4. Three-dimensional Yee cell and electric and magnetic fields staggered in space by half a cell. Electric field components are staggered along their direction, while magnetic field components are staggered perpendicular to their direction (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019).

3. PCE method

Uncertainty Quantification (UQ) quantifies the impact of a system's (the ScaRF code in this work) uncertain input to the system's output quantities of interest (QoI). The PCE method is among the most popular probabilistic UQ methods.

In PCE it is assumed that the ScaRF's input uncertainty is specified by a $d$-dimensional vector $\boldsymbol{\varXi }=(\varXi _{1}, \varXi _{2},\ldots , \varXi _{d} )$ with $\varXi _{i}$, $i=1,\ldots , d$ independent identically distributed variables and realizations $\boldsymbol {\xi }=( \xi _{1}, \xi _{2},\ldots , \xi _{d} )$. Then the joint probability density function (j.p.d.f.) is $f( \boldsymbol {\xi })=\prod _{k=1}^{d}f(\xi _{k} )$, where $f( \xi _{k})$ is the p.d.f. of $\varXi _{k}$. If the QoI $y( \boldsymbol{\varXi} )$ is a scalar with finite variance, $y(\boldsymbol{\varXi})$ can be expanded (PCE method) as

(3.1)\begin{equation} y\left( \boldsymbol{\varXi} \right)=\sum_{j=0}^{\infty} c_{j} \varPsi_{j}\left( \boldsymbol{\varXi}\right),\end{equation}

where $\varPsi _{j}(\boldsymbol{\varXi })$ are multivariate polynomials that are orthogonal with $f( \boldsymbol {\xi } )$ weights, and $c_{j}$ are PCE deterministic coefficients that will be estimated.

In reality the summation in (3.1) is truncated to the first $P+1$ terms:

(3.2)\begin{equation} y\left( \boldsymbol{\varXi} \right)=\sum_{j=0}^{P} c_{j} \varPsi_{j}\left( \boldsymbol{\varXi } \right)+\epsilon \left( \boldsymbol{\varXi} \right),\end{equation}

where $\epsilon ( \boldsymbol{\varXi })$ is the truncation error of the PCE. The construction of the multi-dimensional polynomial basis function $\varPsi _{j}( \boldsymbol{\varXi } )$ is based on the Askey (Xiu & Karniadakis (Reference Xiu and Karniadakis2002)) scheme, which is suitable for general non-Gaussian random inputs. In particular if $\psi _{j_{k}}(\varXi _{k} )$ specify a complete basis set of single variable polynomials of order $j_{k}$, where $j_{k} \in \mathbb {N} \cup \{ 0\}$ and are orthonormal with respect to the p.d.f. $f ( \xi _{k})$, the multi-dimensional basis functions are

(3.3)\begin{equation} \varPsi_{{j}}\left( \boldsymbol{\varXi} \right)=\prod_{k=1}^{d}\psi_{j_{k}}\left( \varXi_{k}\right), \end{equation}

where ${j}=(j_{1},\ldots ,j_{d} )$ specifies the order of $\varPsi _{{j}}( \boldsymbol{\varXi } )$. If the maximum total order of $\varPsi _{{j}}( \boldsymbol{\varXi } )$ is $p_{max}$ ($\sum _{k=1}^{d}j_{k}\leqslant p_{max}$), the number $P$ of basis functions in the truncated PCE expansion (3.2) is

(3.4)\begin{equation} P=\frac{ \left(p_{max}+ d \right)!}{d!p_{max}!}. \end{equation}

If $y$ has finite variance then as $p_{max}$ or $P$ tend to infinity, the PCE in (3.2) is convergent (in the mean-square sense). In addition, if the univariate polynomial basis is chosen based on the Askey scheme, the error term in (3.2) tends to zero exponentially fast (Xiu & Karniadakis Reference Xiu and Karniadakis2002). The PCE coefficients $c_{j}$ are computed by projecting the QoI into the corresponding PCE basis:

(3.5)\begin{equation} c_{j}=E{\left[y\varPsi_{j}\right]}. \end{equation}

In reality, (3.5) is inefficient to apply and so other methods have been developed. These fall into two categories, intrusive and non-intrusive methods. In intrusive methods, the deterministic solvers that relate the system output and input are modified, which can be problematic (Crestaux et al. Reference Crestaux, Le Maitre and Martinez2009). In non-intrusive methods, the deterministic solvers are considered as a black box that provides QoI samples from particular system input samples. In this work, the PCE coefficients are determined using a non-intrusive method based on the Smolyak sparse-grid integration.

3.1. PCE coefficients calculation

The expansion coefficients in (3.2) can be calculated by use of the orthogonality of the basis functions. In particular, the projection of $y$ on each basis function leads to

(3.6)\begin{equation} c_{k} = \frac{\left\langle y,\varPsi_{k} \right\rangle }{\left\| \varPsi_{k} \right\|^{2}},\quad k = 0, \ldots ,P \end{equation}

In (3.6), the numerical calculation of $P+1$ integrals in the numerators need to be performed (the denominators are computed analytically). For this task, the Smolyak sparse-grid integration method is used, because it is more efficient than the full-grid integration techniques, as long as sufficiently smooth quantities are integrated (Crestaux et al. Reference Crestaux, Le Maitre and Martinez2009). The improved performance arises from the fact that sparse grids exhibit similar accuracy levels to full-grid integration with fewer nodes. Therefore, the number of deterministic simulations that are performed is reduced. The function values used in the Smolyak algorithm are located on nodal points of grids given by

(3.7)\begin{equation} \varTheta_{d} = \bigcup_{q=m-d}^{m-1} \,\bigcup_{\boldsymbol{i} \in \mathbb{N}_{q}^{d} }\left( \varTheta_{1}^{i_{1}} \times \varTheta_{1}^{i_{2}} \times \ldots \times \varTheta_{1}^{i_{d}} \right), \end{equation}

where $\varTheta _{1}^{i} = \{\theta _{1}^{i},\theta _{2}^{i},\ldots ,\theta _{m_{i}}^{i} \}$ represent one-dimensional (1-D) sets of $m_{i}$ nodes, $\mathbb {N}_{q}^{d}=\{\boldsymbol {i} \in \mathbb {N}^{d}:\sum _{l=1}^{d}i_{l}=d+q \}$ and $m$ is the accuracy order of the $d$-dimensional quadrature. Actually, the sparse grid is built by combining selected smaller grids that satisfy specific rules. For example, a 2-D sparse grid with accuracy order equal to $3$ is built by considering the $\varTheta _{1}^{2} \times \varTheta _{1}^{1}$, $\varTheta _{1}^{1} \times \varTheta _{1}^{2}$, $\varTheta _{1}^{2} \times \varTheta _{1}^{2}$, $\varTheta _{1}^{3} \times \varTheta _{1}^{1}$ and $\varTheta _{1}^{1} \times \varTheta _{1}^{3}$ grids, and normally includes fewer nodes than the $\varTheta _{1}^{3} \times \varTheta _{1}^{3}$ full grid. After the sparse grid has been defined, the required integrations in (3.6) are calculated by

(3.8)\begin{equation} \left\langle {y,{\varPsi _k}} \right\rangle \simeq \sum_{q=m-d}^{m-1} \sum_{\boldsymbol{i} \in \mathbb{N}_{q}^{d} } {( - 1)}^{m-1-q } \left( {\begin{array}{@{}c@{}}{d - 1} \\ {m-1-q} \end{array}} \right)\left( {{V^{{i_1}}} \otimes {V^{{i_2}}} \otimes \ldots \otimes {V^{{i_d}}}} \right)\left[ y{\varPsi _k} \right] \end{equation}

where the terms $V^{i}[ y\varPsi _{k} ]$ represent the 1-D quadrature. In the case that $y\varPsi _{k}$ is a polynomial of order $2m-1$ or less, then (3.8) is exact.

In this work, the nested Kronrod–Patterson (KP) nodal sets (Petras Reference Petras2003) are used for the construction of the various $\varTheta _1$. These sets are nested, which means that lower-level nodes comprise subsets of higher-level sequences and thus a significant number of re-calculations are avoided when higher level sets are constructed. The number of nodes in a 1-D KP sequence increases at a fast rate, according to

(3.9)\begin{equation} 1, 3, 7, 15, 31, \ldots, 2^{i}-1,\ldots \end{equation}

Therefore, even in the Smolyak algorithm case (3.8), the number of nodes in multidimensional grids grows at an unnecessarily fast rate. To handle this issue, delayed KP sequences are applied (Petras Reference Petras2003) in a way that formulae are reused for some levels before being updated, such that the growth in the number of nodes is under control. In figure 5, the differences between ordinary and delayed KP sequences are shown. In addition, full and sparse 2-D grids are compared in figure 6, where the cases of tensor-based grids, sparse grids with normal KP sequences and sparse grids with delayed KP sequences are examined by considering various construction levels. It is observed that after a few iterations, a tensor grid contains more than $900$ tensor grid nodes. The Smolyak grid with normal KP nodes comprises 129 nodes in the same construction level, while this number reduces to just $33$ nodes for the case where the Smolyak algorithm is applied for delayed sequences.

Figure 5. One-dimensional quadrature nodal distributions corresponding to different accuracy levels. Cases of (a) ordinary and (b) delayed KP sequences are shown (Papadopoulos et al. Reference Papadopoulos, Zygiridis, Glytsis, Kantartzis and Tsiboukis2018).

Figure 6. Comparison of various grids (full grids, normal sparse grids and sparse grids with delayed 1-D sequences) for different accuracy levels (Papadopoulos et al. Reference Papadopoulos, Zygiridis, Glytsis, Kantartzis and Tsiboukis2018).

3.2. Calculation of statistics

The knowledge of the PC expansion coefficients allows the direct computation of the statistical moments of the QoI. The mean value and the standard deviation of $y$ are calculated by

(3.10)\begin{gather} E[y] = c_{0} \end{gather}

(3.11)\begin{gather}\sigma^{2}[y] = \sum_{k = 1}^{P} c_{k}^{2}\left\|\varPsi_{k} \right\|^{2}. \end{gather}

In addition, as will be explained later, sensitivity indices can be easily obtained by the PC coefficients.

3.3. Estimation of the reflection probability density function

Once the PCE coefficients have been calculated, the random output of the system is approximated by the PCE expansion. This easily provides a large enough number of output reflection samples such that a histogram can be built. This is an approximation to the p.d.f. of the reflection. Because a histogram is usually non-smooth, a kernel smoothing method is applied for better visualization of the p.d.f. In particular, the p.d.f. is approximated by

(3.12)\begin{equation} \hat{f}_{R}\left( y \right)=\frac{1}{N_{K}h_{K}}\sum_{i=1}^{N_{K}}K\left( \frac{y-y^{\left( i \right)}}{h_{K}}\right), \end{equation}

where $K(x )$ is a positive-definite function (the kernel) and $h_{K}$ is the bandwidth parameter. A Gaussian kernel (standard normal p.d.f.) and an optimum bandwidth $h_{K}^{\ast }$ (for the Gaussian kernel) is chosen as

(3.13)\begin{equation} h_{K}^{{\ast}}=\left(\frac{2}{3N_{K}}\right)^{1/5} \min\left(\hat{\sigma}_{R}, \hat{iqr}_{R} \right),\end{equation}

where $\hat {\sigma }_{R}$ and $\hat {iqr}_{R}$ are the empirical standard deviation and interquartile, respectively. Here, $\hat {\sigma }_{R}$ is optimal in the sense that it minimizes the $\ell ^{2}$ approximation error $\| \,f_{R}-\hat {f}_{R} \|_{2}$.

3.4. Sensitivity analysis

To quantify the effect of each random input variable $\xi$ (the random amplitude of each mode) the Total Sobol sensitivity indices ($ST_{\xi }$) are calculated (Crestaux et al. Reference Crestaux, Le Maitre and Martinez2009). The $ST_{\xi }$, $\xi \in \{\varLambda ,d_1,d_2,f,\theta \}$ are directly computed from the PC expansion coefficients as

(3.14)\begin{gather} ST_{\xi}=\sum_{u \ni \xi}S_{u} \end{gather}

(3.15)\begin{gather}S_{u}=\frac{\sum_{k \in K_{u}} {c_k^{2}{{\left\| {{\varPsi _k}} \right\|}^{2}}}}{\sum_{k = 1}^{P} {c _k^{2}{{\left\| {{\varPsi _k}} \right\|}^{2}}}} \end{gather}

(3.16)\begin{gather}K_{u}= \left\{k \in \{1,\ldots,P\} \mid {\varPsi _k}(\boldsymbol\xi ) = \prod_{i= 1}^{{\mid} u \mid} {\phi _{\alpha_i^{k}}}({\xi_{u_{k}} } ) , \alpha_{i}^{k}>0.\right\} \end{gather}

where $u \subseteq \{1,2,\ldots ,d \}$. The total Sobol index fpr each random input variable specifies the percentage contribution of the random variable to the variance of the reflection.

4. Numerical results

The ScaRF code, in conjunction with the PCE method described in § 3, is used for the scattering analysis of an incident $O$-mode at 170 GHz (electron cyclotron (EC) frequency range) by a random periodic interface separating blob–background anisotropic media. The interface is a superposition of $5$ spatial sinusoidal modes of wavelengths $c=[0.25, 0.5, 1, 2.5, 5]\lambda _{0}$ and random amplitudes $\xi _{i}$, $i=\{1,\ldots , 5\}$, with a $5\,\%$ variation around $\xi _{i}=1$. The toroidal, $H_{t}$, and poloidal, $H_{p}$, magnetic fields are $4.5T$ and $0.473T$ respectively. The incident wave has a polar angle $\phi =174^{\circ }$ and experiments are performed for samples of azimuthal angles ${\pm }10\,\%$ around $\theta =30^{\circ }$. For the $O$-mode, the normalized wavevector components are $\beta _{x}=0.3496$, $\beta _{y}=0.2007$, $\beta _{z}=-0.0211$, the background and blob densities are $3 \times 10^{20}\ \textrm {m}^{-3}$, $3.2 \times 10^{20}\ \textrm {m}^{-3}$, respectively, and the relative permittivity tensors, (2.1), are calculated via (2.2).

The use of ScaRF allows for the computation of the electric and magnetic field vectors $\boldsymbol {e}$, $\boldsymbol {h}$ at every node of the computational domain, where $\boldsymbol {h}=\textrm {i}Z_{0}^{-1}\boldsymbol {\tilde {h}}$ and therefore the time-averaged Poynting vector, $\boldsymbol {S}$ is computed by

(4.1)\begin{equation} \boldsymbol{S}=\tfrac{1}{2}\mbox{Re}\{\boldsymbol{e} \times \boldsymbol{h}^{*}\}. \end{equation}

Because $\boldsymbol{S}$ is known everywhere in the computational domain, the power reflection $R$, defined as $R=P_{\textrm {ref}}/P_{\textrm {inc}}$, can be calculated simply by numerically integrating $\boldsymbol{S}$ in the $y$ direction, which provides the incident and reflected power $P_{\textrm {inc}}$, $P_{\textrm {ref}}$, respectively. For visualization purposes, in figure 7, the amplitude of the Poynting vector $\boldsymbol{S}$ and the $\boldsymbol{S}$-flow are shown for an incident $O$-mode at $\theta =30^{\circ }$, scattered by a realization of the random density interface.

Figure 7. Amplitude of the Poynting vector ($|\boldsymbol{S}|$) and $\boldsymbol{S}$-flow (red lines: ($S_{x}, S_{y}$) vector) for a plane wave ($O$-mode) incident on a realization of a 5-mode interface (black line) with random amplitudes.

In figure 8, the mean value ($MV$) and standard deviation ($Std$) of $R$ are presented as a function of $\theta$, calculated by the PCE-SG and MC methods. For the PCE-SG and MC methods, 51 and 1000 samples of $R$ were produced by the use of ScaRF, respectively. Therefore the MC method is orders of magnitude more computationally demanding than the PCE-SG method. However, it is observed in figure 8 that the $MV$ and $Std$ values of $R$ are in good agreement, which verifies the accuracy and computational efficiency of the PCE-SG method. The PCE-SG method was used for the maximum polynomial order $p_{\max }=2$ and an accuracy level of the Smolyak grid $m=3$. It is practical and desirable to achieve a minimum value of $R$ so that the maximum power of the incident RF wave will be transmitted into the fusion plasma region. In figure 8, it is concluded that the minimum $MV$ of $R, R\approx 0.04$, occurs at $\theta \approx 30.5^{\circ }$.

Figure 8. Mean Value ($MV$) and Standard Deviation ($Std$) of reflection as a function of the angle of incidence ($\theta$) of the incident wave: (a) calculated by the PCE method; and (b) calculated by the MC method.

The p.d.f. of $R$ is shown in figure 9. It is observed that the maximum of the p.d.f. for $\theta =30.5$ is at $R \approx 0.047$, which, on average, occurs in figure 8 for $\theta \approx 30.5^{\circ }$ as expected.

Figure 9. P.d.f. of reflection for $\theta = 30^{\circ }$.

In figure 10, the Sobol total indices (STI) versus $\theta$ are presented for the five $\xi _{i}$ random variables. The $R$ is more sensitive to the $\xi _{i}$ with the highest STI. From this figure it is concluded that the spatial modes $\lambda _{0}$ and $\lambda _{0}/4$ have respectively the highest and lowest contribution to the variance of $R$, for all angles $\theta$ of the incident RF wave. For the angle $\theta \approx 30.5$, which is of practical interest, the $\lambda _{0}$, $5 \lambda _{0}$ spatial modes have the same contribution to the $R$ variance, the $2.5 \lambda _{0}$, $\lambda _{0} /4$ spatial modes have the lowest contribution to the $R$ variance, while the spatial mode $\lambda _{0} /2$ is somewhere in between.

Figure 10. Sobol Total Indices for the 5 random amplitudes of the spatial modes ($\lambda _{0}/4, \lambda _{0}/2, \lambda _{0}, 2.5\lambda _{0}, 5\lambda _{0}$) that define the density interface.

Finally in figure 11 the 5th and 95th percentiles for the reflection, $R$, are shown. They were calculated by the PCE-SG method, as described in § 3.5, where $N=20\,000$ random reflection samples were used. An interesting observation in this figure is that the desirable minimum reflection $R \approx 0.04$ occurs at $\theta \approx 30.25^{\circ }$ and is a point closer to the $5\,\%$ curve, which means that only roughly $5\,\%$ of the $R$ samples will have reflection smaller than 0.04. Therefore with high probability, launching the incident RF wave at $\theta \approx 30.25^{\circ }$ will result in minimum reflection.

Figure 11. The 5 % and 95 % percentiles calculation as a function of the angle of incidence $\theta$. Percentiles were calculated by the PCE-SG method using $N=20\,000$ reflection samples.

It is noted that the computation of an optimal launch angle of the waves provides guidance to the experimentalists on an appropriate choice of antenna phasing and spectrum. The efficient use of RF waves requires that the coupling of the power to the plasma be optimal, so as to minimize reflections arising from turbulence. Becuase it is not completely possible to quantify the turbulent plasma to any reasonable precision (density variation in space and time), one has to develop statistical tools that can quantify the effect of turbulence and help guide the experimentalists simulations that can be tested on present experiments to gain confidence for future experiments and reactors. The emphasis is not that a certain angle for launching the waves is optimal, but rather there exists such a possibility, that is, there exists a launch window which minimizes the reflection. The existing tools theoretically and/or numerically do not answer the questions: in the presence of a large region of turbulent plasma, what is the most likely scenario for ideal delivery of power and momentum to the core plasma, to minimize reflections and side-scattering, and to control the spectrum of waves propagating in the core so as to optimize the current drive and heating efficiency. One of the first approaches to address these questions is the UQ method of this work whose effectiveness is shown in the presented scattering configuration.