Computation of Electromagnetic Fields Scattered From Objects With Uncertain Shapes Using Multilevel Monte Carlo Method
Alexander Litvinenko, Abdulkadir C. Yucel, Hakan Bagci, Jesper Oppelstrup, Eric Michielssen, Raúl Tempone
CComputation of Electromagnetic Fields Scattered FromOb jects With Uncertain Shapes Using Multilevel MonteCarlo Method
Alexander Litvinenko ∗ , Abdulkadir C. Yucel, Hakan Bagci,Jesper Oppelstrup, Eric Michielssen, and Ra´ul Tempone † September 5, 2018
Abstract
Computational tools for characterizing electromagnetic scattering from objects with un-certain shapes are needed in various applications ranging from remote sensing at microwavefrequencies to Raman spectroscopy at optical frequencies. Often, such computational toolsuse the Monte Carlo (MC) method to sample a parametric space describing geometric uncer-tainties. For each sample, which corresponds to a realization of the geometry, a deterministicelectromagnetic solver computes the scattered fields. However, for an accurate statisticalcharacterization the number of MC samples has to be large. In this work, to address thischallenge, the continuation multilevel Monte Carlo (CMLMC) method is used together witha surface integral equation solver. The CMLMC method optimally balances statistical errorsdue to sampling of the parametric space, and numerical errors due to the discretization ofthe geometry using a hierarchy of discretizations, from coarse to fine. The number of re-alizations of finer discretizations can be kept low, with most samples computed on coarserdiscretizations to minimize computational cost. Consequently, the total execution time issignificantly reduced, in comparison to the standard MC scheme.
Keywords: uncertainty quantification in geometry, random geometry, multilevel Monte Carlomethod (MLMC), continuation MLMC, integral equation, fast multipole method (FMM), fastFourier transform (FFT), scattering problem.
Numerical methods for predicting radar and scattering cross sections (RCS and SCS) of complextargets find engineering applications ranging from microwave remote sensing soil/ocean surface ∗ Corresponding author † A. Litvinenko, H. Bagci, and R. Tempone are with the Strategic Research Initiative - Uncertainty Quan-tification (SRI-UQ) Center, Division of Computer, Electrical, and Mathematical Science and Engineering(CEMSE), King Abdullah University of Science and Technology (KAUST), Thuwal, 23955, Saudi Arabia (e-mails: { alexander.litvinenko, raul.tempone, hakan.bagci } @kaust.edu.sa). A. C. Yucel is with the School of Electricaland Electronics Engineering, Nanyang Technological University, Singapore 639798 (e-mail: [email protected]).J. Oppelstrup is with the Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden(e-mail: [email protected]). E. Michielssen is with the Department of Electrical Engineering and Computer Science,University of Michigan, Ann Arbor, MI 48109, US (e-mail: [email protected]). a r X i v : . [ phy s i c s . c o m p - ph ] S e p nd vegetation [53, 55] to enhancing Raman spectroscopy using metallic nanoparticles [51, 58].When the target size is comparable to or larger than the wavelength at the operation frequency,the scattered field is a strong function of the target shape. However, in many of the applications,whether the target is a (vegetated) rough surface or a nanoparticle, its exact shape may not beknown and has to be parameterized using a stochastic/statistical model. Consequently, compu-tational tools, which are capable of generating statistics of a quantity of interest (QoI) (RCS orSCS in this case) given a geometry described using random variables/parameters, are required.These computational tools often use the Monte Carlo (MC) method [18,27,42,45,48,54,56]. TheMC method is non-intrusive and straightforward to implement; therefore, its incorporation with anexisting deterministic EM solver is rather trivial. However, the traditional MC method has an errorconvergence rate of O ( N − / ) [41], where N is the number of samples in the parametric space usedfor describing the geometry. Provided more regularity of the QoI w.r.t. the geometry parameters,quasi-MC methods may have a better convergence rate, O ( N − ) with a multiplicative log-termthat depends on the number of parameters [29, 41]. Both the traditional MC and quasi-MCmethods require large N to yield accurate statistics of the QoI and are computationally expensivesince the function evaluation at each sampling point, which corresponds to the computation ofthe QoI for one realization of the deterministic problem, requires the execution of a simulation.For practical real-life scattering scenarios, each of these simulations may take a few hours, if nota few days.In recent years, schemes that make use of surrogate models have received significant attention aspotential alternatives to the MC method [14,15,21,31,38,61,62]. The surrogate model is generatedusing the values of the QoI that are computed by the simulator at a small number of “carefullyselected” (collocation) points in the parametric space. The surrogate model is then probed usingthe MC method to obtain statistics of the QoI. Consequently, the surrogate model-based schemesare more efficient than the MC method that operates directly on the computationally expensivesimulator. The efficiency and accuracy of these schemes can be improved using a refinementstrategy that adaptively divides the parametric space into smaller domains, in each of which aseparate collocation scheme is used [14, 15]. Additionally, if the QoI can be approximated interms of low-order contributions from only a few of the parameters, one can use high dimensionalmodel representation expansions to accelerate the generation of the surrogate model [38]. Thesurrogate model-based schemes and and their accelerated variants have been successfully appliedto certain stochastic EM problems [2, 3, 22, 43, 59, 60, 64, 65, 67]. On the other hand, for theproblem of scattering from geometries with uncertain shapes, all parameters and their high-ordercombinations contribute significantly to the QoI limiting the efficient use of a surrogate model [6,7, 69].Recently, the multilevel MC (MLMC) methods have seen increasing use due to their efficiency,robustness, and simplicity [4, 5, 8, 10, 19, 20, 52]. The MLMC methods operate on a hierarchy ofmeshes and perform most of the simulations using low-fidelity models (coarser meshes) and only afew simulations using high-fidelity models (finer meshes). By doing so, their cost becomes signifi-cantly smaller than the cost of the traditional MC methods using only high-fidelity models. Thecontinuation MLMC (CMLMC) algorithm [10] uses the MLMC method to calibrate the averagecost per sampling point and the corresponding variance using Bayesian estimation, taking partic-ular notice of the deepest levels of the mesh hierarchy, to minimize the computational cost. Tobalance discretization and statistical errors, the CMLMC method estimates how the discretizationerror and the computational cost depend on the mesh level and uses this information to selectoptimal numbers of levels and samples on each level.2 x y x´ z´ y´ Illumination: (0,0,-1) Rotation: ( ϕ x , ϕ y , ϕ z ) Observation: ( θ,ϕ ) Figure 1: Description of the scattering problem.In this work, a computational framework, which makes use of the CMLMC method to ef-ficiently and accurately characterize EM wave scattering from dielectric objects with uncertainshapes, is proposed. The deterministic simulations required by the CMLMC method to computethe samples at different levels are carried out using the Poggio-Miller-Chan-Harrington-Wu-Tsaisurface integral equation (PMCHWT-SIE) solver [39]. The PMCHWT-SIE is discretized usingthe method of moments (MoM) and the iterative solution of the resulting matrix system is ac-celerated using a (parallelized) fast multipole method (FMM) - fast Fourier transform (FFT)scheme (FMM-FFT) [9, 11, 49, 50, 57, 66, 68]. Numerical results, which demonstrate the accuracy,efficiency, and convergence of the proposed computational framework, are presented. The devel-oped computational framework is proven effective not only in studying the effects of uncertaintyin the geometry on scattered fields but also in increasing the robustness of the FMM-FFT accel-erated PMCHWT-SIE solver by testing its convergence for a large set of scenarios with deformedgeometries and varying mesh densities, quadrature rules, iterative solver tolerances, and FMMparameters. 3
Formulation
This section describes a computational framework for characterizing scattering from dielectric ob-jects with uncertain shapes (1). It is assumed that the shape of the scatterer can be parameterizedby means of random variables. The QoI is the SCS over a user-defined solid angle (i.e., a measureof far-field scattered power in a cone). Section 2.1 describes the MLMC and CMLMC methods,Section 2.2 formulates a scheme to parameterize and generate objects with uncertain shapes, andSection 2.3 outlines the FMM-FFT accelerated PMCHWT-SIE solver used for computing thescattering cross section of a given object.
This section describes elements of the MLMC and CMLMC methods relevant to the character-ization of scattering from objects with uncertain shapes. A more in-depth description of thesetechniques can be found in [10].Let ξ and g ( ξ ) represent the vector of random parameters/variables and the QoI, respectively.The goal of the MLMC method is to approximate the expected value, E[ g ], to a guaranteedtolerance TOL with predefined probability, and minimal computational cost. To achieve this,the method constructs a telescoping sum, defined over a sequence of mesh levels (cid:96) = 0 , . . . , L asdescribed next.Let { P (cid:96) } L(cid:96) =0 and { h (cid:96) } L(cid:96) =0 be sequences of meshes that discretize the randomly generated object’ssurface and their average element edge sizes, respectively. It is assumed { P (cid:96) } L(cid:96) =0 are generatedhierarchically with h (cid:96) = h β − (cid:96) for h > β >
1. A method to such meshes for thepurpose of illustrating the proposed technique is described in Section 2.2.Let g (cid:96) ( ξ ) represent the approximation to g ( ξ ) computed using mesh P (cid:96) . The MLMC methodexpresses E[ g L ], the expected value of the most accurate approximation g L , usingE[ g L ] = L (cid:88) (cid:96) =0 E[ G (cid:96) ] (1)where G (cid:96) is defined as G (cid:96) = (cid:40) g if (cid:96) = 0 g (cid:96) − g (cid:96) − if (cid:96) > . (2)Note that g (cid:96) and g (cid:96) − are computed using the same input random parameter ξ .In the telescoping sum (1), the expected values in practice are replaced by sample averages,i.e., E[ G (cid:96) ] ≈ ∼ G (cid:96) = M − (cid:96) (cid:80) M (cid:96) m =1 G (cid:96),m , where random variable G (cid:96),m have the same distribution as G (cid:96) and are independent identically distributed (i.i.d.) samples. As (cid:96) increases, the variance of G (cid:96) decreases. As a result, the total computational cost can be reduced by approximating E[ G (cid:96) ] witha smaller number of samples.The CMLMC algorithm is an improved version of the MLMC method in that it approximatesE[ g ] with a sequence of decreasing tolerances [10]. In doing so, it continuously improves estimatesof several problem-dependent parameters, while solving relatively inexpensive problems that bythemselves would yield large tolerances. The CMLMC algorithm assumes that the convergence4ates for the mean (weak convergence) and variance (strong convergence) followE[ g − g (cid:96) ] ≈ Q W h q (cid:96) (3a)Var[ g (cid:96) − g (cid:96) − ] ≈ Q S h q (cid:96) − (3b)for Q W (cid:54) = 0, Q S > q >
0, and 0 < q ≤ q [10]. For example, the CMLMC algorithm estimates q ≈ q ≈
4, and q ≈ q ≈ ‘cost persample’ and ‘total cost’ as described next.The CMLMC estimator for the QoI, A , can be written as A = (cid:80) L(cid:96) =0 ∼ G (cid:96) . Let the average costof generating one sample of G (cid:96) (cost of one deterministic simulation for one random realization)be W (cid:96) ∝ h − dγ(cid:96) = h − dγ β (cid:96)dγ (4)where d is the spatial dimension and γ is determined by the computational complexity of thedeterministic solver. For the FMM-FFT-accelerated PMCHWT-SIE solver used here d = 2 (onlysurfaces are discretized) and γ ≈ G (cid:96) may fluctuate for different realizations depending on the number of iterationsrequired. Finally, for the method used for generating { P (cid:96) } Ll =0 (Section 2.2), β = 2.The total CMLMC computational cost is W = L (cid:88) (cid:96) =0 M (cid:96) W (cid:96) . (5)The estimator A satisfies a tolerance with a prescribed failure probability 0 < ν ≤
1, i.e.,P[ | E[ g ] − A| ≤ TOL] ≥ − ν (6)while minimizing W . The total error is split into bias and statistical error, | E[ g ] − A| ≤ | E[ g − A ] | (cid:124) (cid:123)(cid:122) (cid:125) Bias + | E[ A ] − A| (cid:124) (cid:123)(cid:122) (cid:125) Statistical error where θ ∈ (0 ,
1) is a splitting parameter, so thatTOL = (1 − θ )TOL (cid:124) (cid:123)(cid:122) (cid:125) Bias tolerance + θ TOL (cid:124) (cid:123)(cid:122) (cid:125)
Statistical error tolerance . (7)The CMLMC algorithm bounds the bias, B = | E[ g − A ] | , and the statistical error as B = | E[ g − A ] | ≤ (1 − θ )TOL (8) | E[ A ] − A| ≤ θ TOL (9)5here the latter bound holds with probability 1 − ν . Note that θ itself can be a variable [10].To satisfy condition in (9) we require:Var[ A ] ≤ (cid:18) θ TOL C ν (cid:19) (10)for some given confidence parameter, C ν , such that Φ( C ν ) = 1 − ν , (see more in [23]); here, Φ isthe cumulative distribution function of a standard normal random variable.By construction of the MLMC estimator, E[ A ] = E[ g L ], and by independence, Var[ A ] = (cid:80) L(cid:96) =0 V (cid:96) M − (cid:96) , where V (cid:96) = Var[ G (cid:96) ]. Given L , TOL, and 0 < θ <
1, and by minimizing W subjectto the statistical constraint (10) for { M (cid:96) } L(cid:96) =0 gives the following optimal number of samples perlevel (cid:96) (apply ceiling function to M (cid:96) if necessary): M (cid:96) = (cid:18) C ν θ TOL (cid:19) (cid:114) V (cid:96) W (cid:96) (cid:32) L (cid:88) (cid:96) =0 (cid:112) V (cid:96) W (cid:96) (cid:33) . (11)Summing the optimal numbers of samples over all levels yields the following expression for thetotal optimal computational cost in terms of TOL: W (TOL , L ) = (cid:18) C ν θ TOL (cid:19) (cid:32) L (cid:88) (cid:96) =0 (cid:112) V (cid:96) W (cid:96) (cid:33) . (12)The total cost of the CMLMC algorithm can be estimated using Theorem 1 below [5, 8, 19, 24, 25]. Theorem 1
Let d = { , , } denote the problem dimension. Suppose there exist positive con-stants q , q , γ > such that q ≥ min( q , γd ) , and | E[ g (cid:96) − g ] | = O ( h q (cid:96) ) , Var[ g (cid:96) − g (cid:96) − ] = O ( h q (cid:96) ) ,and W (cid:96) = O ( h − dγ(cid:96) ) . Then for any accuracy TOL and confidence level ν , < ν ≤ , there exist adeepest level L ( TOL ) and number of realizations { M (cid:96) ( TOL ) } such that lim TOL → inf P[( | E ( g ) − A| ≤ TOL )] ≥ − ν (13) with cost W ( TOL ) ≤ O (cid:0) TOL − (cid:1) , q > dγ O (cid:16) TOL − (cid:0) log(TOL − ) (cid:1) (cid:17) , q = dγ O (cid:18) TOL − (cid:16) dγ − q q (cid:17) (cid:19) , q < dγ . (14)This theorem shows that even in the worst case scenario, the CMLMC algorithm has alower computational cost than that of the traditional (single level) MC method, which scalesas O (TOL − − dγ/q ). Furthermore, in the best case scenario presented above, the computationalcost of the CMLMC algorithm scales as O (cid:0) TOL − (cid:1) , i.e. identical to that of the MC methodassuming that the cost per sample is O (1) . In other words, for this case, the CMLMC algorithmcan in effect remove the computational cost required by the discretization, namely O (TOL − dγ/q ).6 .2 Random Shape Generation This section describes a method to generate the sequence of meshes { P (cid:96) } L(cid:96) =0 that are used inthe numerical experiments (Section 3) to demonstrate the properties of the CMLMC scheme.Alternative ways to generating random perturbations can be found in [28, 30, 32–34]. The methodused here relies on perturbing an initial mesh generated on a sphere and refining it as level (cid:96) isincreased as described next.First, the unit sphere is discretized using a mesh P with triangular elements, then the pertur-bations are generated by moving the nodes of this mesh using v ( ϑ m , ϕ m ) ≈ ˜ v ( ϑ m , ϕ m ) + K (cid:88) k =1 a k κ k ( ϑ m , ϕ m ) . (15)where ϑ m and ϕ m are angular coordinates of node m , v ( ϑ m , ϕ m ) is its (perturbed) radial coordinate,and ˜ v ( ϑ m , ϕ m ) = 1 m is its (unperturbed) radial coordinate on the unit sphere (all in sphericalcoordinate system). Here, κ k ( ϑ, ϕ ) are obtained from spherical harmonics by re-scaling theirarguments and a k , which satisfy (cid:80) Kk =1 a k < .
5, are uncorrelated random variables.For the numerical experiments considered in this work, K = 2, and κ ( ϑ, ϕ ) = cos( α ϑ ) and κ ( ϑ, ϕ ) = sin( α ϑ ) sin( α ϕ ), where α , α , and α are positive constants. If a κ k ( ϑ, ϕ ) dependson ϕ , then its dependence on ϑ must vanish at the poles ϑ = { , π } . Therefore, α must be aninteger. Additionally, α must be an integer to generate a smooth perturbation.The random surface generated using (15) is also rotated and scaled as described next. Let¯ R x ( ϕ x ), ¯ R y ( ϕ y ), and ¯ R z ( ϕ z ) represent the matrices defined to perform rotations around axes x , y , and z by angles ϕ x , ϕ y , and ϕ z , respectively:¯ R x ( ϕ x ) = ϕ x − sin ϕ x ϕ x cos ϕ x ¯ R y ( ϕ y ) = cos ϕ y ϕ y − sin ϕ y ϕ y ¯ R z ( ϕ z ) = cos ϕ z − sin ϕ z ϕ z cos ϕ z
00 0 1 . Similarly, let ¯ L ( l x , l y , l z ) represent the matrix defined to implement scaling along axes x , y , and z by l x , l y , and l z , respectively: ¯ L ( l x , l y , l z ) = /l x /l y
00 0 1 /l z . Then, the coordinates for the “rotated” and “scaled” mesh P are obtained after the applicationof (15), rotation and scaling matrices: x (cid:48) m y (cid:48) m z (cid:48) m = ¯ L ( l x , l y , l z ) ¯ R x ( ϕ x ) ¯ R y ( ϕ y ) ¯ R z ( ϕ z ) x m y m z m . (16)7ere, ( x m , y m , z m ), and ( x (cid:48) m , y (cid:48) m , z (cid:48) m ) are the coordinates of node m before and after the transfor-mation.The random variables used in generating the final version of P are the perturbation weights a k , k = 1 , . . . , K , the rotation angles ϕ x , ϕ y , and ϕ z , and the scaling factors l x , l y , and l z , makingthe dimension of the stochastic space K + 6, i.e., random parameter vector ξ = { a , . . . , a K , ϕ x , ϕ y , ϕ z , l x , l y , l z } . Note that P is the coarsest mesh used in CMLMC ( (cid:96) = 0) (see Fig. 2(a) for an example).The mesh of the next level ( (cid:96) = 1), P , is generated by refining each triangle of the perturbed P into four (by halving all three edges and connecting mid-points) [Fig. 2(b)]. The mesh at level (cid:96) = 2, P , is generated in the same way from P [Fig. 2(c)], and so on. All meshes P (cid:96) at all levels (cid:96) = 1 , . . . , L are nested discretizations of P . This method of refinement results in β = 2 in (4).Note that no uncertainties are added on meshes P (cid:96) , (cid:96) >
0; the uncertainty is introduced only atlevel (cid:96) = 0. It is assumed that P is fine enough to accurately represent the variations of thehighest harmonic in the expansion used in (15). This section briefly describes the PMCHWT-SIE solver used for the RCS and SCS of a dielectricscatterer with surface S . It is assumed that sources and fields are time-harmonic, i.e., their timedependence varies as e jωt , where t and ω are the time and angular frequency, respectively.Let V and V represent the space internal and external to S , respectively. The permittivity,permeability, and characteristic impedance in V i , i ∈ { , } are ε i , µ i , and η i = (cid:112) µ i /ε i , re-spectively. It is assumed that the scatterer is excited by an external electromagnetic field withelectric and magnetic components E inc ( r ) and H inc ( r ). Using the surface equivalence theoremand enforcing the tangential continuity of total electromagnetic fields on S yield the PMCHWTSIE [39]: ˆn ( r ) × E inc ( r ) = ˆn ( r ) × (cid:110) L [ J ]( r ) + L [ J ]( r ) (17) − K [ M ]( r ) − K [ M ]( r ) (cid:111) , r ∈ S ˆn ( r ) × H inc ( r ) = ˆn ( r ) × (cid:110) K [ J ]( r ) + K [ J ]( r ) (18)+ η − L [ M ]( r ) + η − L [ M ]( r ) (cid:111) , r ∈ S. Here, ˆn ( r ) is the outward pointing unit normal vector, and J ( r ) and M ( r ) represent equivalentelectric and magnetic surface current densities on S . The integral operators L i [ . ]( r ) and K i [ . ]( r )in (17) and (18) are L i [ X ]( r ) = jωµ i (cid:90) S (cid:20) ¯I + ∇∇ (cid:48) k i (cid:21) · X ( r (cid:48) ) g i ( r , r (cid:48) ) ds (cid:48) K i [ X ]( r ) = ∇ × (cid:90) S X ( r (cid:48) ) g i ( r , r (cid:48) ) ds (cid:48) a) (b)(c) (d) Figure 2: An example of four nested meshes with (a) 320, (b) 1280, (c) 5120, and (d) 20480triangular elements, which are generated for a perturbed shape with α = 2, α = 3, α = 2, a = 0 .
04 m, a = 0 .
048 m, ϕ x = 0 .
32 rad, ϕ y = 0 .
88 rad, ϕ z = 0 .
81 rad, l x = 1 . l y = 1 .
08, and l z = 1 . g i ( r , r (cid:48) ) = e − jk i | r − r (cid:48) | / (4 π | r − r (cid:48) | ) is the Green function of the Helmholtz equation in theunbounded medium with wave number k i = ω √ ε i µ i .To numerically solve (17) and (18) for the unknowns J ( r ) and M ( r ), S is discretized by9riangular elements, and J ( r ) and M ( r ) are approximated as J ( r ) = N (cid:88) n =1 ¯ I n f n ( r ) (19) M ( r ) = N (cid:88) n = N +1 ¯ I n f n − N ( r ) (20)where f n ( r ), n = 1 , ..., N , represent Rao-Wilton-Glisson basis functions [44], and I = [ I , . . . , I N ] T is the vector of unknown coefficients. Substituting (19) and (20) into (17) and (18), and testingthe resulting equations with f m ( r ), m = 1 , ..., N , yield the MoM system of equations:¯ V = ¯ Z ¯ I (21)where the entries of the vector ¯ V and the MOM matrix ¯ Z are¯ V m = (cid:40)(cid:10) f m , E inc (cid:11) , ≤ m ≤ N (cid:10) f m − N , H inc (cid:11) , N + 1 ≤ m ≤ N (22)¯ Z m,n = (cid:104) f m , L [ f n ] + L [ f n ] (cid:105) , ≤ m ≤ N , ≤ n ≤ N − (cid:104) f m − N , K [ f n ] + K [ f n ] (cid:105) ,N + 1 ≤ m ≤ N , ≤ n ≤ N (cid:104) f m , K [ f n − N ] + K [ f n − N ] (cid:105) , ≤ n ≤ N , N + 1 ≤ n ≤ N (cid:10) f m − N , η − L [ f n ] + η − L [ f n ] (cid:11) ,N + 1 ≤ m ≤ N , N + 1 ≤ n ≤ N . (23)Here, the inner product is (cid:104) f m , a (cid:105) = (cid:82) S f m ( r ) · a ( r ) d r .Matrix equation (21) is solved iteratively for ¯ I . The computational cost of multiplying ¯ Z witha trial solution vector scales as O ( N iter N ), where N iter is the number of iterations required forthe residual error to reach the desired level: typically N iter (cid:28) N . Likewise, the storage costs ofthe unaccelerated/classical solution scale as O ( N ).To minimize the computational and storage cost while executing the MLMC algorithm a fastmultipole method (FMM)-fast Fourier transform (FFT) scheme is used. A detailed formulation ofFMM and its extension FMM-FTT can be found in [9,11,49,50,57,66,68]. The scheme encloses thescatterer in a fictitious box that is embedded into a uniform grid of smaller boxes. Two non-emptyboxes (i.e., boxes containing at least a pair of patches) constitute a far-field pair if there is at leastone box between them. Otherwise, they form a near-field pair. Interactions between basis and testfunctions in near-field pairs (and same boxes) are computed using (23) and stored in matrix ¯ Z near .The contribution of near- and self-interactions to the matrix-vector multiplication ¯ Z ¯ I is calculatedby simply multiplying ¯ Z near with ¯ I . The interactions between basis and test functions in far-fieldpairs are represented using the radiation patterns of basis and test functions (sampled over asolid angle of a sphere) and a translation operator and their contributions to ¯ Z ¯ I are computedusing the FMM-FFT scheme. The computational cost and memory requirement of the FMM-FFT10ccelerated iterative solution scale as O ( N iter N / log / N ) and O ( N / log / N ), respectively [66,68]. Note that these estimates are higher than those of the multilevel FMM [46,47]; the FMM-FFTscheme is preferred here since its parallel implementation is significantly simpler [12, 13, 16, 35–37,40, 49, 50, 57, 66, 68]. The parallelization strategy implemented here uses a hybrid message passinginterface/open multiprocessing (MPI/OpenMP) standard to uniformly distribute the memory andcomputational load among processors [12, 13, 16, 35–37, 40, 49, 50, 57, 66, 68].Note that the computational cost estimate provided above is obtained under the assumptionthat the ratio of the wavelength to the (average) element edge length stays the same as N isdecreased/increased [49, 50, 57, 66, 68]. Within the CMLMC algorithm, the mesh is refined fromone level to the next while the frequency is kept constant. This means that the ratio of thewavelength to the edge length increases and the computational cost of the FFT-FMM-acceleratedsolver is expected to scale differently. Indeed, numerical experiments in Section 3 show that theCMLMC algorithm estimates the computational cost parameter as γ ≈ E inc ( r ) = E e − jk ˆu inc · r and H inc ( r ) = ( ˆu inc × E ) /η e − jk ˆu inc · r , where ˆu inc is the direction of propagation. The unknownvector ¯ I is solved for under this excitation and the RCS σ rcs ( ϑ, ϕ ) along direction ˆu ( ϑ, ϕ ) = ˆx sin ϑ cos ϕ + ˆy sin ϑ sin ϕ + ˆz cos ϑ is computed using [26] σ rcs ( ϑ, ϕ ) = (cid:12)(cid:12) F ( ϑ, ϕ ) (cid:12)(cid:12) π (cid:12)(cid:12) E (cid:12)(cid:12) . (24)Here, F ( ϑ, ϕ ) is the scattered electric field pattern in the far field and computed using F ( ϑ, ϕ ) = − jωµ ˆu ( ϑ, ϕ ) × ˆu ( ϑ, ϕ ) × N (cid:88) n =1 ¯ I n (cid:90) S n f n ( r (cid:48) ) e jk ˆu ( ϑ,ϕ ) · r (cid:48) d r (cid:48) − jk ˆu ( ϑ, ϕ ) × N (cid:88) n = N +1 ¯ I n (cid:90) S n − N f n − N ( r (cid:48) ) e jk ˆu ( ϑ,ϕ ) · r (cid:48) d r (cid:48) . The SCS C sca (Ω) is obtained by integrating σ rcs ( ϑ, ϕ ) over the solid angle Ω [26]: C sca (Ω) = 14 π (cid:90) Ω σ rcs ( ϑ, ϕ ) sin ϑdϑdϕ. (25)The integral in (25) is efficiently computed using the exact quadrature rule in [63]. In all numerical experiments considered in this section, the scatterer resides in free space (vacuum)with µ = 4 π × − H/m and ε = c /µ F/m, where c is the speed of light in vacuum, andthe scatterer’s permittivity and permeability are µ = µ and ε = 4 ε , respectively. For theplane wave excitation, E = ˆx , ˆu inc = − ˆ z , and ω = 2 πf with f = 300 MHz. The QoI is thescattering cross-section C sca (Ω) computed over the cone Ω = [ ϑ , ϑ ] × [ ϕ , ϕ ] = [1 / , / π rad × [5 / , / π rad. For all random perturbations, the parameters of the quasi-spherical harmonics11 a) (b)(c) (d) Figure 3: Amplitudes of (a) J ( r ) and (b) M ( r ) induced on the unit sphere under excitation byan ˆ x -polarized plane wave propagating in − ˆ z direction at 300 MHz. Amplitudes of (c) J ( r ) and(d) M ( r ) induced on the perturbed shape shown in Fig. 2(c) under excitation by the same planewave. For all figures, amplitudes are normalized to 1 and plotted in dB scale.12 /4 /2 /4 0 /4 /2 3 /4 (rad) -10-50510152025 r c s ( d B ) SpherePerturbed surface (a) (rad) -10-50510152025 r c s ( d B ) SpherePerturbed surface (b)
Figure 4: RCS of the unit sphere and the perturbed shape shown in Fig. 2(c) computed on (a) xz and (b) yz planes under excitation by an ˆ x -polarized plane wave propagating in − ˆ z directionat 300 MHz. For (a) ϕ = 0 and ϕ = π rad in the first and second halves of the horizontalaxis, respectively. For (b), ϕ = π/ ϕ = 3 π/ α = 2, α = 3, and α = 2. Meshes P (cid:96) at levels (cid:96) = 0 , , , ,
4, which are generated usingthe method described in Section 2.2, have { , , , , } triangles, respectively.Note that not all five mesh levels are required in every experiment. Since the CMLMC algorithmis stochastic, 15 independent CMLMC runs are executed for each experiment to statisticallycharacterize the performance of the method.The matrix system in (21) is solved iteratively using the transpose-free quasi-minimal residual(TFQMR) method [17] with a tolerance (residual) of (cid:15) (cid:96) = (cid:15) β − (cid:96) , where (cid:96) = 0 , , , , (cid:15) = 6 . × − , and β = 2 (Section 2.2). Using a level-dependent tolerance isconsistent with the dependence of the discretization error in the QoI, which scales at a (heuristicallydetermined) rate O ( h (cid:96) ) ∝ O (cid:0) β − (cid:96) (cid:1) and ensures that the number of iterations is kept in checkwhen analyzing coarser meshes. The integrals in (22) and (23) are computed using the Gaussianquadrature rules; the accuracy of the rule, i.e. number of quadrature points, adjusts with the meshlevel: it is seven for (cid:96) = 0 , , (cid:96) = 3 ,
4. For all levels of meshes, the parameters of theFMM-FFT scheme are selected carefully to ensure that it has six digits of accuracy [46, 47]. Allsimulations are executed on Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80 GHz DELL workstationswith 40 cores and 128 GB RAM. The electromagnetic simulator is implemented in Fortran 90while the CMLMC is written in C++ and Python [1]. The interface between the electromagneticsolver and MLMC library is implemented in C.
This section demonstrates that the scattered fields strongly depend on the shape of the object forthe scenarios considered in the numerical experiments, i.e., f = 300 MHz, ε = 4 ε , and the sizeof the object is roughly 2 m. More specifically, J ( r ) and M ( r ) induced on a (rotated and scaled)perturbed geometry and its RCS are compared to the same quantities of the unit sphere. Therotated, scaled, and perturbed surface is generated using a = 0 .
04 m, a = 0 .
048 m, ϕ x = 0 .
32 rad, ϕ y = 0 .
88 rad, ϕ z = 0 .
81 rad, l x = 1 . l y = 1 .
08, and l z = 1 .
07. The discretization is refinedtwice, resulting in the mesh with 5120 triangles [Fig. 2(c)]. A mesh with the same number oftriangles is generated on the unit sphere. Figs 3(a) and 3(b) show normalized amplitudes of J ( r )13nd M ( r ) induced on this unit sphere. Fig. 3(c) and 3(d) plot the normalized amplitudes of J ( r )and M ( r ) on the perturbed shape. It is clear that there is a significant difference between thecurrent distributions induced on the sphere and the perturbed shaped. Figs. 4(a) and 4(b) compare σ rcs ( ϑ, ϕ ) of the unit sphere and the perturbed shaped computed on the xz − ( ϑ ∈ [0 , π ] rad, ϕ = 0and ϕ = π rad) and yz − ( ϑ ∈ [0 , π ] rad, ϕ = π/ ϕ = 3 π/ θ . Additionally, performances of the CMLMCand MC methods are compared with theoretical estimates. Since the number of mesh points, N ,is changing from level to level, the computational time and computational work, is shown in termsof TOL, rather than in N (e.g., O (TOL − )). In this example, the CMLMC algorithm is executed for random variables a , a ∼ U [ − . , . ϕ x , ϕ y , ϕ z ∼ U [0 . ,
3] rad, and l x , l y , l z ∼ U [0 . , . U [ a, b ] is the uniform distributionbetween a and b . The CMLMC algorithm is run for TOL values ranging from 0 . . P (cid:96) , (cid:96) = 0 , , , , − − − TOL A v e r a g e T i m e ( s ) TOL − CMLMCMC Estimate (a) − − − TOL W o r k e s t i m a t e TOL − CMLMCMC Estimate (b)
Figure 5: (a) Computation times required by the CMLMC and MC methods vs. TOL. (b) Valueof the computational cost estimate W [given by (5)] for the CMLMC and MC methods vs. TOL.Computation times are averaged over 15 repetitions of the experiment.Fig. 5(a) compares computation times for the CMLMC and MC algorithms as a function ofTOL. The figure reveals three convergence regimes (TOL zones).The first zone, 0 . (cid:47) TOL (cid:47) .
3, covers the range where TOL is larger than the sum of bias andstatistical error. No additional samples and no additional mesh refinements are required by theCMLMC algorithm. The MC and CMLMC methods consume similar computational resources.14he second zone, 0 . (cid:47) TOL (cid:47) .
1, corresponds to a pre-asymptotic regime. Both con-vergence rates are approximately 2. The statistical error is dominant, and many new samples onexisting meshes P (cid:96) , (cid:96) = 0 , , (cid:96) ).The third zone, 0 . (cid:47) TOL (cid:47) . (cid:39) . P (cid:96) , (cid:96) = 0 , , , (cid:47) .
012 anadditional finer mesh, P (cid:96) , (cid:96) = 4, is required. Sampling on level (cid:96) = 4 is more expensive; therefore,the MC computation time increases rapidly for TOL (cid:47) . (cid:96) = 4, and is, therefore, significantly faster.The computational cost estimate W is an indicator of computation time. It depends on howthe computational cost of the deterministic solver changes from level (cid:96) − (cid:96) [as indicated byparameters γ and β in (4)] and on the order of decay for the mean and the variance [parameters q , q in (3a)]. Fig. 5(b) plots the values of W vs. TOL. The curve is similar to that of theCMLMC computation time given in Fig. 5(a) demonstrating that γ ≈ q ≈
3, and q ≈ G (cid:96) = g (cid:96) − g (cid:96) − vs. (cid:96) . Computation times varyroughly as 2 (cid:96) , which verifies that γ ≈ d = 2 and β = 2 in (4)]. Fig. 6(b) shows E (cid:96) = E[ G (cid:96) ]vs. (cid:96) , revealing that E (cid:96) ∼ − (cid:96) (assumed weak convergence obtained with q = 3). Fig. 6(c) shows V (cid:96) = Var[ G (cid:96) ] vs (cid:96) , demonstrating that V (cid:96) ∼ − (cid:96) (assumed strong convergence with q = 5).The results presented in Figs. 5(a) and 5(b), and Figs. 6(a), 6(b), and 6(c) confirm theassumptions stated in Section 2.1 as well as the CMLMC scheme’s quasi-optimality.15 ‘ T i m e ( s ) ‘ G ‘ (a) ‘ − − − − − − E ‘ − ‘ G ‘ (b) ‘ − − − − − − − V ‘ − ‘ G ‘ (c) Figure 6: (a) Time required to compute G (cid:96) vs. (cid:96) . (b) E (cid:96) = E[ G (cid:96) ] vs. (cid:96) and assumed weakconvergence curve 2 − (cid:96) ( q = 3). (c) V (cid:96) = Var[ G (cid:96) ] vs. (cid:96) and assumed strong convergence curve2 − (cid:96) ( q = 5). The experiment is repeated 15 times independently and the obtained values areshown as error bars on the curves. 16 − − − TOL . . . . . . θ Figure 7: Value of θ used by the CMLMC algorithm vs. TOL. Computation time is averagedover 15 repetitions of the experiment, i.e., there are 15 values of θ for a given value of TOL.Fig. 7 shows θ vs. TOL and demonstrates the complex relationship between TOL and thebias and statistical error. θ decreases from 1 to ≈ . θ = 1 implies that the ratio of bias tothe total error is negligible, i.e. that the ratio of of statistical error to the total error is 1. Suchvariability in θ is one of the differences between the CMLMC algorithm and the MLMC methodwhere θ = 0 .
5. The blue dots show that for TOL ≈ . θ ≈
1, meaning that the impact of thebias is negligible and that there is no need to further extend the mesh hierarchy by adding finermesh levels.A higher accuracy can be achieved by decreasing either the bias or the statistical error (i.e.,smaller TOL). θ ≈ ≈ . G (cid:96) is to look at its probability density function (pdf).Figs. 8(a) and 8(b) plot empirical pdfs of G (cid:96) from { , , } samples for (cid:96) = 1 and (cid:96) = { , } ,respectively. Fig. 8(a) shows that Var[ G ] = Var[ g − g ], where g and g are computed usingthe meshes P and P , varies roughly in the range (0 , . g − g ] ≈ .
0. Fig. 8(b) shows that G = g − g , where g are g are computed using themeshes P and P , varies in the interval ( − . , .
1) and E[ g − g ] ≈ .
02. Finally, G = g − g ,where g are g are computed using the meshes P and P , varies in the interval ( − . , . G ] ≈
0. The pdfs of G (cid:96) concentrate more and more around zero [with the rates shown inFigs. 6(b), and 6(c)] as (cid:96) increases. In this example, the CMLMC algorithm is executed for random variables a , a ∼ U [0 . , . ϕ x , ϕ y , ϕ z ∼ U [0 . , .
0] rad, and l x , l y , l z ∼ U [0 . , . g1-g0 (a) -0.2 -0.1 0 0.1 0.2020406080100120 g2-g1g3-g2 (b) Figure 8: Pdfs of g (cid:96) − g (cid:96) − for (a) (cid:96) = 1 and (b) (cid:96) = { , } .changed from 0 . .
02. For the lowest TOL value, the CMLMC algorithm requires four meshlevels, i.e., P (cid:96) , (cid:96) = 0 , , , . (cid:47) TOL (cid:47) .
15, describes the regimewhen TOL is higher than the sum of the bias and statistical error. No additional samples orrefinements are needed. − − TOL A v e r a g e T i m e ( s ) TOL − CMLMCMC Estimate (a) − − TOL W o r k e s t i m a t e TOL − CMLMCMC Estimate (b)
Figure 9: (a) Computation times required by the CMLMC and MC methods vs. TOL. (b) Valueof the computational cost estimate W [given by (4)] for the CMLMC and MC methods vs. TOL.Computations times are averaged over 15 repetitions of the experiment.The second zone, TOL (cid:47) .
07, describes the pre-asymptotic regime. As TOL gets smaller,the CMLMC algorithm becomes more efficient than the MC method. For values of TOL closeto 0 .
02, the CMLMC algorithm is roughly 10 times faster than MC. Fig. 9(b) shows the valuesof the computational cost estimate W vs. TOL. The curve is similar to that of the CMLMCcomputation time given in Fig. 9(a) validating the model developed for the CMLMC algorithm.18igs. 10(a) and 10(b) show E (cid:96) = E[ G (cid:96) ] and V (cid:96) = Var[ G (cid:96) ] vs. (cid:96) , respectively, revealing de-pendencies on (cid:96) that vary as 2 − (cid:96) (assuming weak convergence obtained with q = 2) and 2 − (cid:96) (assuming strong convergence obtained with q = 4). Note that in the previous example theseparameters are 3 and 5, respectively, demonstrating that the dependence of q and q on thevariability of the random variables used to described the geometry. ‘ − − − − E ‘ − ‘ G ‘ (a) ‘ − − − − − − − V ‘ − ‘ G ‘ (b) Figure 10: (a) E (cid:96) = E[ G (cid:96) ] vs. (cid:96) and assumed weak convergence curve 2 − (cid:96) ( q = 2). (b) V (cid:96) = Var[ G (cid:96) ] vs. (cid:96) and assumed strong convergence curve 2 − (cid:96) ( q = 4). The experiment isrepeated 15 times and error bars are shown on the curves.19 Conclusion
A computational framework is developed to efficiently and accurately characterize EM wave scat-tering from dielectric objects with uncertain shapes. To this end, the framework uses the CMLMCalgorithm, which reduces the computational cost of the traditional MC method by performing mostof the simulations with lower accuracy and lower cost (using coarser meshes) and smaller numberof simulations with higher accuracy and higher cost (using finer meshes). To increase the efficiencyfurther, each of the simulations is carried out using the FMM-FFT accelerated PMCHWT-SIEsolver. Numerical results demonstrate that the CMLMC algorithm can be 10 times faster thanthe traditional MC method depending on the amplitude of the perturbations used for representingthe uncertainties in the scatterer’s shape. This work confirms that the known advantages of theCMLMC algorithm can be observed when it is applied to EM wave scattering: non-intrusiveness,dimension independence, better convergence rates compared to the classical MC method, andhigher immunity to irregularity w.r.t. uncertain parameters, than, for example, sparse grid meth-ods.For optimal performance (for a given value of accuracy parameter TOL), the CMLMC algo-rithm requires the mean and the variance to have reliable convergence rates (i.e., one should beable estimate q and q without much difference from one level to next). However, some randomperturbations may affect the convergence rates. With difficult-to-predict convergence rates, it ishard for the CMLMC algorithm to estimate the computational cost W , the number of levels L ,the number of samples on each level M (cid:96) , the computation time, and the parameter θ , and thevariance in QoI. All these may result in a sub-optimal performance. Indeed, numerical resultsdemonstrate that there is a significant pre-asymptotic regime where the performance is not opti-mal. Additionally, it is observed that the settings of the FMM-FFT accelerated PMCHWT-SIEsolver, which regulate the computation time and the accuracy (such as the iterative solver thresh-old, the number of quadrature points, and the FMM-FFT parameters), have a significant effecton the performance of the CMLMC algorithm.It should also be noted here that the CMLMC algorithm is proven effective not only in studyingthe effects of uncertainty in the geometry on scattered fields but also in increasing the robustnessof the FMM-FFT accelerated PMCHWT-SIE solver by testing its convergence for a large set ofscenarios with deformed geometries and varying mesh densities, quadrature rules, iterative solvertolerances, and FMM parameters. Acknowledgements
The research reported in this publication was supported by funding from King Abdullah Universityof Science and Technology (KAUST). We gratefully acknowledge support from Abdul-Lateef Haji-Ali (KAUST and Oxford), developer of the CMLMC, Erik von Schwerin and Haakon Hoel forvaluable comments, concerning the CMLMC algorithm, as well as Ismail Enes Uysal and H¨useyinArda ¨Ulk¨u (KAUST) for their initial efforts on this project.
References [1] https://github.com/StochasticNumerics/mimclib.git .202] A. C. M. Austin, N. Sood, J. Siu, and C. D. Sarris. Application of polynomial chaos to quantifyuncertainty in deterministic channel models.
IEEE Trans. Antennas Propag. , 61(11):5754–5761, Nov. 2013.[3] H. Bagci, A. C. Yucel, J. S. Hesthaven, and E. Michielssen. A fast Stroud-based collocationmethod for statistically characterizing EMI/EMC phenomena on complex platforms.
IEEETrans. Electromagn. Compat. , 51(2):301–311, May 2009.[4] A. Barth, C. Schwab, and N. Zollinger. Multi-level Monte Carlo finite element method forelliptic PDEs with stochastic coefficients.
Numer. Math. , 119(1):123–161, 2011.[5] J. Charrier, R. Scheichl, and A. L. Teckentrup. Finite element error analysis of elliptic pdeswith random coefficients and its application to multilevel Monte Carlo methods.
SIAM J.Numer. Anal. , 51(1):322–352, 2013.[6] C. Chauvi´ere, J. S. Hesthaven, and L. Lurati. Computational modeling of uncertainty intime-domain electromagnetics.
SIAM J. Sci. Comput. , 28(2):751–775, 2006.[7] C. Chauvi´ere, J. S. Hesthaven, and L. C. Wilcox. Efficient computation of RCS from scatterersof uncertain shapes.
IEEE Trans. Electromagn. Compat. , 55(5):1437–1448, May 2007.[8] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup. Multilevel Monte Carlo methodsand applications to elliptic PDEs with random coefficients.
Comput. Vis. Sci. , 14(1):3–15,2011.[9] R. Coifman, V. Rokhlin, and S. Wandzura. The fast multipole method for the wave equation:a pedestrian prescription.
IEEE Antennas Propag. Mag. , 35(3):7–12, June 1993.[10] N. Collier, A.-L. Haji-Ali, F. Nobile, E. von Schwerin, and R. Tempone. A continuationmultilevel Monte Carlo algorithm.
BIT Numer. Math. , 55(2):399–432, 2015.[11] N. Engheta, W. D. Murphy, V. Rokhlin, and M. S. Vassiliou. The fast multipole method(FMM) for electromagnetic scattering problems.
IEEE Trans. Antennas Propag. , 40(6):634–641, June 1992.[12] O. Ergul and L. Gurel. Efficient parallelization of the multilevel fast multipole algorithm forthe solution of large-scale scattering problems.
IEEE Trans. Antennas Propag. , 56(8):2335–2345, Aug. 2008.[13] O. Ergul and L. Gurel. A hierarchical partitioning strategy for an efficient parallelizationof the multilevel fast multipole algorithm.
IEEE Trans. Antennas Propag. , 57(6):1740–1750,June 2009.[14] J. Foo and G. E. Karniadakis. Multi-element probabilistic collocation method in high dimen-sions.
J. Comput. Phys. , 229(5):1536–1557, 2010.[15] J. Foo, X. Wan, and G. E. Karniadakis. The multi-element probabilistic collocation method(ME-PCM): Error analysis and applications.
J. Comput. Phys. , 227(22):9572–9595, 2008.[16] J. Fostier and F. Olyslager. An asynchronous parallel MLFMA for scattering at multipledielectric objects.
IEEE Trans. Antennas Propag. , 56(8):2346–2355, Aug. 2008.2117] R. W. Freund. A transpose-free quasi-minimal residual algorithm for non-Hermitian linearsystems.
SIAM J. Sci. Stat. Comput. , 14:470–482, Mar. 1993.[18] N. Garcia and E. Stoll. Monte Carlo calculation for electromagnetic-wave scattering fromrandom rough surfaces.
Phys. Rev. Lett. , 52:1798–1801, May 1984.[19] M. B. Giles. Multilevel Monte Carlo path simulation.
Operations Research , 56(3):607–617,2008.[20] M. B. Giles. Multilevel Monte Carlo methods.
Acta Numerica , 24:259–328, 2015.[21] L. Giraldi, A. Litvinenko, D. Liu, H. G. Matthies, and A. Nouy. To be or not to be intrusive?the solution of parametric and stochastic equations—the “plain vanilla” Galerkin case.
SIAMJ. Sci. Comput. , 36(6):A2720–A2744, 2014.[22] L. J. Gomez, A. C. Ycel, L. Hernandez-Garcia, S. F. Taylor, and E. Michielssen. Uncertaintyquantification in transcranial magnetic stimulation via high-dimensional model representa-tion.
IEEE Trans. Biomed. Eng. , 62(1):361–372, Jan. 2015.[23] A.-L. Haji-Ali, F. Nobile, E. von Schwerin, and R. Tempone. Optimization of mesh hierarchiesin multilevel Monte Carlo samplers.
Stoch. Partial Differ. Equ. Anal. Comput. , 4(1):76–112,Mar. 2016.[24] H. Hoel, E. Von Schwerin, A. Szepessy, and R. Tempone. Adaptive multilevel Monte Carlosimulation. In
Numerical Analysis of Multiscale Computations , pages 217–234. Springer, 2012.[25] H. Hoel, E. Von Schwerin, A. Szepessy, and R. Tempone. Implementation and analysisof an adaptive multilevel Monte Carlo algorithm.
Monte Carlo Methods and Applications ,20(1):1–41, 2014.[26] A. Ishimaru.
Wave Propagation and Scattering in Random Media . Academic Press, NewYork, US, 1978.[27] V. Jandhyala, E. Michielssen, S. Balasubramaniam, and W. C. Chew. A combined steepestdescent-fast multipole algorithm for the fast analysis of three-dimensional scattering by roughsurfaces.
IEEE Trans. Geosci. Remote Sens. , 36(3):738–748, May 1998.[28] B. N. Khoromskij, A. Litvinenko, and H. G. Matthies. Application of hierarchical matricesfor computing the Karhunen-Lo`eve expansion.
Computing , 84(1-2):49–67, 2009.[29] F. Kuo, R. Scheichl, C. Schwab, I. Sloan, and E. Ullmann. Multilevel quasi-monte carlomethods for lognormal diffusion problems.
Mathematics of Computation , 86(308):2827–2860,2017.[30] A. Lang and C. Schwab. Isotropic Gaussian random fields on the sphere: Regularity, fastsimulation and stochastic partial differential equations.
ArXiv e-prints , May 2013.[31] A. Litvinenko, H. Matthies, and T. A. El-Moselhy. Sampling and low-rank tensor approxi-mation of the response surface. In J. Dick, F. Y. Kuo, G. W. Peters, and I. H. Sloan, editors,
Monte Carlo and Quasi-Monte Carlo Methods 2012 , volume 65 of
Springer Proceedings inMathematics & Statistics , pages 535–551. Springer Berlin Heidelberg, 2013.2232] A. Litvinenko and H. G. Matthies. Sparse data representation of random fields.
Proc. Appl.Math. Mech. , 9(1):587–588, 2009.[33] A. Litvinenko and H. G. Matthies. Numerical methods for uncertainty quantification andbayesian update in aerodynamics. In
Management and Minimisation of Uncertainties andErrors in Numerical Aerodynamics , pages 265–282. Springer, 2013.[34] D. Liu, A. Litvinenko, C. Schillings, and V. Schulz. Quantification of airfoil geometry-inducedaerodynamic uncertainties—comparison of approaches.
SIAM/ASA J. Uncertain. Quantif. ,5(1):334–352, 2017.[35] Y. Liu, A. Al-Jarro, H. Bagci, and E. Michielssen. Parallel PWTD-accelerated explicit so-lution of the time-domain electric field volume integral equation.
IEEE Trans. AntennasPropag. , 64(6):2378–2388, June 2016.[36] Y. Liu, A. C. Yucel, H. Bagci, A. C. Gilbert, and E. Michielssen. A wavelet-enhancedPWTD-accelerated time-domain integral equation solver for analysis of transient scatteringfrom electrically large conducting objects.
IEEE Trans. Antennas Propag. , 66(5):2458–2470,May 2018.[37] Y. Liu, A. C. Yucel, H. Bagci, and E. Michielssen. A scalable parallel PWTD-accelerated siesolver for analyzing transient scattering from electrically large objects.
IEEE Trans. AntennasPropag. , 64(2):663–674, Feb. 2016.[38] X. Ma and N. Zabaras. An adaptive high-dimensional stochastic model representationtechnique for the solution of stochastic partial differential equations.
J. Comput. Phys. ,229(10):3884–3915, 2010.[39] L. N. Medgyesi-Mitchang, J. M. Putnam, and M. B. Gedera. Generalized method of momentsfor three-dimensional penetrable scatterers.
J. Opt. Soc. Am. , A11:1383–1398, 1994.[40] B. Michiels, J. Fostier, I. Bogaert, and D. D. Zutter. Weak scalability analysis of thedistributed-memory parallel MLFMA.
IEEE Trans. Antennas Propag. , 61(11):5567–5574,Nov. 2013.[41] W. J. Morokoff and R. E. Caflisch. Quasi-Monte Carlo integration.
J. Comput. Phys. ,122(2):218–230, 1995.[42] M. Nieto-Vesperinas and J. M. Soto-Crespo. Monte Carlo simulations for scattering of electro-magnetic waves from perfectly conductive random rough surfaces.
Opt. Express , 12(12):979–981, July 1987.[43] J. S. Ochoa and A. C. Cangellaris. Random-space dimensionality reduction for expedientyield estimation of passive microwave structures.
IEEE Trans. Microwave Theory Tech. ,61(12):4313–4321, Dec. 2013.[44] S. Rao, D. Wilton, and A. Glisson. Electromagnetic scattering by surfaces of arbitrary shape.
IEEE Trans. Antennas Propag. , 30(3):409–418, May 1982.[45] K. Sarabandi, P. F. Polatin, and F. T. Ulaby. Monte Carlo simulation of scattering from alayer of vertical cylinders.
IEEE Trans. Antennas Propag. , 41(4):465–475, Apr. 1993.2346] J. Song and W. C. Chew. Multilevel fast multipole algorithm for solving combined fieldintegral equations of electromagnetic scattering.
Microw. Opt. Technol. Lett. , 10(1):14–19,1995.[47] J. Song, C.-C. Lu, and W. C. Chew. Multilevel fast multipole algorithm for electromagneticscattering by large complex objects.
IEEE Trans. Antennas Propag. , 45(10):1488–1493, Oct.1997.[48] J. M. Soto-Crespo and M. Nieto-Vesperinas. Electromagnetic scattering from very roughrandom surfaces and deep reflection gratings.
J. Opt. Soc. Am. A , 6(3):367–384, Mar. 1989.[49] J. M. Taboada, M. G. Araujo, F. O. Basteiro, J. L. Rodriguez, and L. Landesa. MLFMA-FFT parallel algorithm for the solution of extremely large problems in electromagnetics.
Proc.IEEE , 101:350–363, 2013.[50] J. M. Taboada, L. Landesa, F. Obelleiro, J. L. Rodriguez, J. M. Bertolo, M. G. Araujo, J. C.Mourino, and A. Gomez. High scalability FMM-FFT electromagnetic solver for supercom-puter systems.
IEEE Antennas Propag. Mag. , 51:20–28, 2009.[51] C. E. Talley, J. B. Jackson, C. Oubre, N. K. Grady, C. W. Hollars, S. M. Lane, T. R. Huser,P. Nordlander, and N. J. Halas. Surface-enhanced Raman scattering from individual Aunanoparticles and nanoparticle dimer substrates.
Nano Lett. , 5(8):1569–1574, 2005.[52] A. L. Teckentrup, R. Scheichl, M. B. Giles, and E. Ullmann. Further analysis of multilevelMonte Carlo methods for elliptic PDEs with random coefficients.
Numer. Math. , 125(3):569–600, Nov. 2013.[53] L. Tsang, K. Ding, S. Huang, and X. Xu. Electromagnetic computation in scattering ofelectromagnetic waves by random rough surface and dense media in microwave remote sensingof land surfaces.
Proc. IEEE , 101(2):255–279, Feb. 2013.[54] L. Tsang, K. H. Ding, S. E. Shih, and J. A. Kong. Scattering of electromagnetic waves fromdense distributions of spheroidal particles based on Monte Carlo simulations.
J. Opt. Soc.Am. A , 15(10):2660–2669, Oct. 1998.[55] L. Tsang, T. Liao, S. Tan, H. Huang, T. Qiao, and K. Ding. Rough surface and volumescattering of soil surfaces, ocean surfaces, snow, and vegetation based on numerical Maxwellmodel of 3-D simulations.
IEEE J. Sel. Topics Appl. Earth Observ. , 10(11):4703–4720, Nov.2017.[56] R. L. Wagner, J. Song, and W. C. Chew. Monte Carlo simulation of electromagnetic scatteringfrom two-dimensional random rough surfaces.
IEEE Trans. Antennas Propag. , 45(2):235–245,Feb. 1997.[57] C. Waltz, K. Sertel, M. A. Carr, B. C. Usner, and J. L. Volakis. Massively parallel fast mul-tipole method solutions of large electromagnetic scattering problems.
IEEE Trans. AntennasPropag. , 55(6):1810–1816, June 2007.[58] H. Wang, C. S. Levin, and N. J. Halas. Nanosphere arrays with controlled sub-10-nm gaps assurface-enhanced Raman spectroscopy substrates.
J. Am. Chem. Soc. , 127(43):14992–14993,2005. 2459] T.-W. Weng, Z. Zhang, Z. Su, Y. Marzouk, A. Melloni, and L. Daniel. Uncertainty quan-tification of silicon photonic devices with correlated and non-Gaussian random parameters.
Opt. Lett. , 23(4):4242–4254, Feb. 2015.[60] Y. Xing, D. Spina, A. Li, T. Dhaene, and W. Bogaerts. Stochastic collocation for device-levelvariability analysis in integrated photonics.
Photon Res. , 4(2):93–100, Apr. 2016.[61] D. Xiu. Efficient collocational approach for parametric uncertainty analysis.
Commun. Com-put. Phys. , 2(2):293–309, Apr. 2007.[62] D. Xiu and J. S. Hesthaven. High-order collocation methods for differential equations withrandom inputs.
SIAM J. Sci. Comput. , 27(3):1118–1139, 2005.[63] A. C. Yucel. Helmholtz and high-frequency Maxwell multilevel fast multipole algorithms withself-tuning library. Master’s thesis, University of Michigan, Ann Arbor, Michigan, US, 2008.[64] A. C. Yucel, H. Bagci, and E. Michielssen. An adaptive multi-element probabilistic colloca-tion method for statistical EMC/EMI characterization.
IEEE Trans. Electromagn. Compat. ,55(6):1154–1168, Dec. 2013.[65] A. C. Yucel, H. Bagci, and E. Michielssen. An ME-PC enhanced HDMR method for efficientstatistical analysis of multiconductor transmission line networks.
IEEE Trans. Compon.Packag. Manuf. Technol. , 5(5):685–696, May 2015.[66] A. C. Yucel, L. J. Gomez, and E. Michielssen. Compression of translation operator tensorsin FMM-FFT-accelerated SIE solvers via Tucker decomposition.
IEEE Antennas WirelessPropag. Lett. , 16:2667–2670, 2017.[67] A. C. Yucel, Y. Liu, H. Bagci, and E. Michielssen. Statistical characterization of electro-magnetic wave propagation in mine environments.
IEEE Antennas Wireless Propag. Lett. ,12:1602–1605, 2013.[68] A. C. Yucel, W. Sheng, C. Zhou, Y. Liu, H. Bagci, and E. Michielssen. An FMM-FFTaccelerated SIE simulator for analyzing EM wave propagation in mine environments loadedwith conductors.
IEEE J. Multiscale and Multiphysics Comput. Techn. , 3:3–15, 2018.[69] Z. Zeng and J.-M. Jin. Efficient calculation of scattering variation due to uncertain geometricaldeviation.