Download as PDF, PPTX
![CMLMC numerical tests The QoI is the SCS over a user-defined solid angle Ω = [1/6, 11/36]π rad × [5/12, 19/36]π rad (i.e., a measure of far-field scattered power in a cone). RVs are: a1, a2 ∼ U[−0.14, 0.14] m, ϕx, ϕy, ϕz ∼ U[0.2, 3] rad, lx, ly, lz ∼ U[0.9, 1.1]; U[a, b] is the uniform distribution between a and b. CMLMC runs for TOL ranging from 0.2 to 0.008. At TOL ≈ 0.008, CMLMC requires L = 5 meshes with {320, 1280, 5120, 20480, 81920} triangles. 15/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-16-2048.jpg)
![V = Var[G ] vs. (strong convergence) 0 1 2 3 4 10−7 10−6 10−5 10−4 10−3 10−2 10−1 V 2−5 G assumed strong convergence curve 2−5 (q2 = 5). 20/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-21-2048.jpg)
![CMLMC repeating Let {P }L =0 be sequences of meshes with h = h0β− , β > 1. Let g (ξ) represent the approximation to g(ξ) computed using mesh P . E[gL] = L =0 E[G ] (4) where G is defined as G = g0 if = 0 g − g −1 if > 0 . (5) Note that g and g −1 are computed using the same input random parameter ξ. 28/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-29-2048.jpg)
![CMLMC repeating E[G ] ≈ ∼ G = M−1 M m=1 G ,m, E[g − g ] ≈ QW hq1 (6a) Var[g − g −1] ≈ QShq2 −1 (6b) for QW = 0, QS > 0, q1 > 0, and 0 < q2 ≤ 2q1. QoI A = L =0 ∼ G . Let the average cost of generating one sample of G (cost of one deterministic simulation for one random realization) be W ∝ h−dγ = h−dγ 0 β dγ (7) 29/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-30-2048.jpg)
![CMLMC repeating The total CMLMC computational cost is W = L =0 M W . (8) The estimator A satisfies a tolerance with a prescribed failure probability 0 < ν ≤ 1, i.e., P[|E[g] − A| ≤ TOL] ≥ 1 − ν (9) while minimizing W . The total error is split into bias and statistical error, |E[g] − A| ≤ |E[g − A]| Bias + |E[A] − A| Statistical error 30/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-31-2048.jpg)
![CMLMC repeating Let θ ∈ (0, 1) be a splitting parameter, so that TOL = (1 − θ)TOL Bias tolerance + θTOL Statistical error tolerance . (10) The CMLMC algorithm bounds the bias, B = |E[g − A]|, and the statistical error as B = |E[g − A]| ≤ (1 − θ)TOL (11) |E[A] − A| ≤ θTOL (12) where the latter bound holds with probability 1 − ν. To satisfy condition in (12) we require: Var[A] ≤ θTOL Cν 2 (13) for some given confidence parameter, Cν, such that Φ(Cν) = 1 − ν 2, Φ is the cdf of a standard normal random variable. 31/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-32-2048.jpg)
![CMLMC repeating By construction of the MLMC estimator, E[A] = E[gL], and by independence, Var[A] = L =0 V M−1 , where V = Var[G ]. Given L, TOL, and 0 < θ < 1, and by minimizing W subject to the statistical constraint (13) for {M }L =0 gives the following optimal number of samples per level (apply ceiling function to M if necessary): M = Cν θTOL 2 V W L =0 V W . (14) Summing the optimal numbers of samples over all levels yields the following expression for the total optimal computational cost in terms of TOL: W (TOL, L) = Cν θTOL 2 L =0 V W 2 . (15) 32/33](https://image.slidesharecdn.com/computationofelectromagneticfieldsscatteredfromdielectricobjectsofuncertainshapesusingmlmc-200306133636/75/Computation-of-electromagnetic-fields-scattered-from-dielectric-objects-of-uncertain-shapes-using-MLMC-algorithm-33-2048.jpg)
Computation of electromagnetic fields scattered from dielectric objects of uncertain shapes using MLMC algorithm
- 1. Computation of ElectromagneticFields Scattered From Dielectric Objects of Uncertain Shapes Using MLMC A. Litvinenko1 , A. C. Yucel3 , H. Bagci2 , J. Oppelstrup4 , E. Michielssen5 , R. Tempone1,2 1 RWTH Aachen, 2 KAUST, 3 Nanyang Technological University in Singapore, 4 KTH Royal Institute of Technology, 5 University of Michigan GAMM 2020 in Kassel, Germany
- 2. Motivation: Comput. toolsfor electromagnetic scattering from objects with uncertain shapes are needed in various applications. Goal: Develop numerical methods for predicting radar and scattering cross sections (RCS and SCS) of complex targets. How: To reduce comp. cost we use CMLMC method (advanced version of Multi-Level Monte Carlo). CMLMC optimally balances statistical and discretization errors. It requires very few samples on fine meshes and more on coarse. taken from wiki, reddit.com, EMCoS 1/33
- 3. Plan: 1. Scattering problemsetup 2. Deterministic solver 3. Generation of random shapes 4. Shape transformation 5. QoI on perturbed shape 6. Continuation Multi Level Monte Carlo (CMLMC) 7. Results (time, work vs. TOL, weak and strong convergences) 8. Conclusion 9. Best practices 10. Appendix (details of CMLMC) 2/33
- 4. Scattering problem Excitation: Observation: x y z (r),0 0 ,0 (,) ˆu(,) ˆuinc (inc ,inc ) (inc ,inc ) Input:randomly perturbed shape Output: radar and scattering cross sections, electric and magnetic surface current densities 3/33
- 5. Deterministic solver Electromagnetic scatteringfrom dielectric objects is analyzed by using the Poggio-Miller-Chan-Harrington-Wu-Tsai surface integral equation (PMCHWT-SIE) solver. The PMCHWT-SIE is discretized using the method of moments (MoM) and the iterative solution of the resulting matrix system is accelerated using a (parallelized) fast multipole method (FMM) - fast Fourier transform (FFT) scheme (FMM-FFT). 4/33
- 6. Generation of randomshapes Perturbed shape v(ϑm, ϕm) is defined as v(ϑm, ϕm) ≈ ˜v(ϑm, ϕm) + K k=1 akκk(ϑm, ϕm). (1) where ϑm and ϕm are angular coordinates of node m, ˜v(ϑm, ϕm) = 1 m is unperturbed radial coordinate on the unit sphere. κk(ϑ, ϕ) obtained from spherical harmonics by re-scaling their arguments, κ1(ϑ, ϕ) = cos(α1ϑ), κ2(ϑ, ϕ) = sin(α2ϑ) sin(α3ϕ), where α1, α2, α3 > 0. -1 1 -0.5 1 0 0.5 0 1 0 -1 -1 1 0 -1-1 -0.5 0 0.5 1 0.5 0 -0.5 -1 1 1 0 -1 1 0 -1 1 0.5 0 -0.5 -1 -1.5 5/33
- 7. Mesh transformation Mesh P0,which is now after the application of (1), is also rotated and scaled using the simple transformation xm ym zm :=L(lx, ly, lz)Rx(ϕx)Ry(ϕy)Rz(ϕz) xm ym zm . (2) Node (xm, ym, zm) before and after transformation. matrices Rx(ϕx), Ry(ϕy), and Rz(ϕz) perform rotations around x, y, and z axes by angles ϕx, ϕy, and ϕz, matrix L(lx, ly, lz) implements scaling along x, y, and z axes by lx, ly, and lz, respectively. 6/33
- 8. Random rotation, stretchingand expanding rotations around axes x, y, and z by angles ϕx, ϕy, and ϕz: Rx(ϕx) = 1 0 0 0 cos ϕx − sin ϕx 0 sin ϕx cos ϕx Ry(ϕy) = cos ϕy 0 sin ϕy 0 1 0 − sin ϕy 0 cos ϕy Rz(ϕz) = cos ϕz − sin ϕz 0 sin ϕz cos ϕz 0 0 0 1 . Matrix L(lx, ly, lz) implements scaling along axes x, y, z by factors lx, ly, and lz: ¯L(lx, ly, lz) = 1/lx 0 0 0 1/ly 0 0 0 1/lz . 7/33
- 9. Input random vector RVsused in generating the coarsest perturbed mesh P0 are: 1. perturbation weights ak, k = 1, . . . , K, 2. rotation angles ϕx, ϕy, and ϕz, 3. scaling factors lx, ly, and lz. Thus, random input parameter vector: ξ = (a1, . . . , aK , ϕx, ϕy, ϕz, lx, ly, lz) ∈ RK+6 (3) defines the perturbed shape. 8/33
- 10. Mesh refinement Mesh P1( = 1) is generated by refining each triangle of the perturbed P0 into four (by halving all three edges and connecting mid-points). Mesh at level = 2, P2, is generated in the same way from P1. All meshes P at all levels = 1, . . . , L are nested discretizations of P0. (!!!) No uncertainties are added on meshes P , > 0; the uncertainty is introduced only at level = 0. 9/33
- 11. Refinement of theperturbed shape 4 nested meshes with {320, 1280, 5120, 20480} triangular elements. 10/33
- 12. Electric (left) andmagnetic (right) surface current densities 11/33
- 13. Electric (left) andmagnetic (right) surface current densities Amplitudes of (a) J(r) and (b) M(r) induced on the unit sphere under excitation by an ˆ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 under excitation by the same plane wave. For all figures, amplitudes are normalized to 1 and plotted in dB scale. 12/33
- 14. RCS of unitsphere and perturbed shape 3 /4 /2 /4 0 /4 /2 3 /4 (rad) -10 -5 0 5 10 15 20 25 rcs (dB) Sphere Perturbed surface 3 /4 /2 /4 0 /4 /2 3 /4 (rad) -10 -5 0 5 10 15 20 25 rcs (dB) Sphere Perturbed surface RCS is computed on (a) xz and (b) yz planes under excitation by an ˆx-polarized plane wave propagating in −ˆz di- rection at 300 MHz. (a) ϕ = 0 and ϕ = π rad in the first and second halves of the hori- zontal axis, respectively. (b)ϕ = π/2 rad and ϕ = 3π/2 rad in the first and second halves of the horizontal axis. 13/33
- 15. Multilevel Monte CarloAlgorithm Aim: to approximate the mean E (g(u)) of QoI g(u) to a given accuracy ε := TOL, where u = u(ω) and ω - random perturbations in the domain. Idea: Balance discretization and statistical errors. Very few samples on a fine grid and more on coarse (denote by M ). Assume: have a hierarchy of L + 1 meshes {h }L =0, h := h0β− for each realization of random domain. 14/33
- 16. CMLMC numerical tests TheQoI is the SCS over a user-defined solid angle Ω = [1/6, 11/36]π rad × [5/12, 19/36]π rad (i.e., a measure of far-field scattered power in a cone). RVs are: a1, a2 ∼ U[−0.14, 0.14] m, ϕx, ϕy, ϕz ∼ U[0.2, 3] rad, lx, ly, lz ∼ U[0.9, 1.1]; U[a, b] is the uniform distribution between a and b. CMLMC runs for TOL ranging from 0.2 to 0.008. At TOL ≈ 0.008, CMLMC requires L = 5 meshes with {320, 1280, 5120, 20480, 81920} triangles. 15/33
- 17. Average time vs.TOL 10−3 10−2 10−1 100 TOL 104 105 106 107 108 AverageTime(s) TOL−2 CMLMC MC Estimate 16/33
- 18. Work estimate vs.TOL 10−3 10−2 10−1 100 TOL 101 102 103 104 105 106 Workestimate TOL−2 CMLMC MC Estimate 17/33
- 19. Time required tocompute G vs. . 0 1 2 3 4 101 102 103 104 105 Time(s) 22 G 18/33
- 20. E = E(G ) vs. (weak convergence) 0 1 2 3 4 10−6 10−5 10−4 10−3 10−2 10−1 100 E 2−3 G assumed weak convergence curve 2−3 (q1 = 3). 19/33
- 21. V = Var[G] vs. (strong convergence) 0 1 2 3 4 10−7 10−6 10−5 10−4 10−3 10−2 10−1 V 2−5 G assumed strong convergence curve 2−5 (q2 = 5). 20/33
- 22. Value of θ 10−3 10−2 10−1 100 TOL 0.0 0.2 0.4 0.6 0.8 1.0 θ Theexperiment is repeated 15 times independently and the obtained values are shown as error bars on the curves. 21/33
- 23. Prob. density functionsof g − g −1 -1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 g1-g0 -0.2 -0.1 0 0.1 0.2 0 20 40 60 80 100 120 g2-g1 g3-g2 (a) = 1 and (b) = {2, 3}. 22/33
- 24. Best practices forapplying CMLMC method to CEM problems Download CMLMC: https://github.com/StochasticNumerics/mimclib.git (or use MLMC from M. Giles) Implement interface to couple CMLMC and your deterministic solver Generate a hierarchy of meshes (mimimum 3), nested are better Generate 5-7 random shapes on first 3 meshes Estimate the strong and weak convergence rates, q1, q2, (later they will be corrected by CMLMC algorithm) Run CMLMC solver and check visually automatically generated plots 23/33
- 25. Conclusion Used CMLMC methodto characterize EM wave scattering from dielectric objects with uncertain shapes. Researched how uncertainties in the shape propagate to the solution. Demonstrated that the CMLMC algorithm can be 10 times faster than MC. To increase the efficiency further, each of the simulations is carried out using the FMM-FFT accelerated PMCHWT-SIE solver. Confirmed that the known advantages of the CMLMC 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, and higher immunity to irregularity w.r.t. uncertain parameters, than, for example, sparse grid methods. 24/33
- 26. Conclusion WARING: Some randomperturbations may affect the convergence rates in CMLMC. With difficult-to-predict convergence rates, it is hard for the CMLMC algorithm to estimate the computational cost W , the number of levels L, the number of samples on each level M , the computation time, and the parameter θ, and the variance in QoI. All these may result in a sub-optimal performance. 25/33
- 27. Acknowledgements SRI UQ atKAUST and Alexander von Humboldt foundation. Results are published: 1. A. Litvinenko, A. C. Yucel, H. Bagci, J. Oppelstrup, E. Michielssen, R. Tempone, Computation of Electromagnetic Fields Scattered From Objects With Uncertain Shapes Using Multilevel Monte Carlo Method, IEEE J. on Multiscale and Multiphysics Comput. Techniques, pp 37-50, 2019. 2. A. Litvinenko, A. C. Yucel, H. Bagci, J. Oppelstrup, E. Michielssen, R. Tempone, Computation of Electromagnetic Fields Scattered From Objects With Uncertain Shapes Using Multilevel Monte Carlo Method, arXiv:1809.00362, 2018 26/33
- 28. Appendix Let us repeatthe main idea of (C)MLMC method 27/33
- 29. CMLMC repeating Let {P}L =0 be sequences of meshes with h = h0β− , β > 1. Let g (ξ) represent the approximation to g(ξ) computed using mesh P . E[gL] = L =0 E[G ] (4) where G is defined as G = g0 if = 0 g − g −1 if > 0 . (5) Note that g and g −1 are computed using the same input random parameter ξ. 28/33
- 30. CMLMC repeating E[G ]≈ ∼ G = M−1 M m=1 G ,m, E[g − g ] ≈ QW hq1 (6a) Var[g − g −1] ≈ QShq2 −1 (6b) for QW = 0, QS > 0, q1 > 0, and 0 < q2 ≤ 2q1. QoI A = L =0 ∼ G . Let the average cost of generating one sample of G (cost of one deterministic simulation for one random realization) be W ∝ h−dγ = h−dγ 0 β dγ (7) 29/33
- 31. CMLMC repeating The totalCMLMC computational cost is W = L =0 M W . (8) The estimator A satisfies a tolerance with a prescribed failure probability 0 < ν ≤ 1, i.e., P[|E[g] − A| ≤ TOL] ≥ 1 − ν (9) while minimizing W . The total error is split into bias and statistical error, |E[g] − A| ≤ |E[g − A]| Bias + |E[A] − A| Statistical error 30/33
- 32. CMLMC repeating Let θ∈ (0, 1) be a splitting parameter, so that TOL = (1 − θ)TOL Bias tolerance + θTOL Statistical error tolerance . (10) The CMLMC algorithm bounds the bias, B = |E[g − A]|, and the statistical error as B = |E[g − A]| ≤ (1 − θ)TOL (11) |E[A] − A| ≤ θTOL (12) where the latter bound holds with probability 1 − ν. To satisfy condition in (12) we require: Var[A] ≤ θTOL Cν 2 (13) for some given confidence parameter, Cν, such that Φ(Cν) = 1 − ν 2, Φ is the cdf of a standard normal random variable. 31/33
- 33. CMLMC repeating By constructionof the MLMC estimator, E[A] = E[gL], and by independence, Var[A] = L =0 V M−1 , where V = Var[G ]. Given L, TOL, and 0 < θ < 1, and by minimizing W subject to the statistical constraint (13) for {M }L =0 gives the following optimal number of samples per level (apply ceiling function to M if necessary): M = Cν θTOL 2 V W L =0 V W . (14) Summing the optimal numbers of samples over all levels yields the following expression for the total optimal computational cost in terms of TOL: W (TOL, L) = Cν θTOL 2 L =0 V W 2 . (15) 32/33
- 34. Literature 1. Collier, N.,Haji-Ali, A., Nobile, F. et al. A continuation multilevel Monte Carlo algorithm. Bit Numer Math 55, 399–432 (2015). https://doi.org/10.1007/s10543-014-0511-3 2. D Liu, A Litvinenko, C Schillings, V Schulz, Quantification of Airfoil Geometry-Induced Aerodynamic Uncertainties—Comparison of Approaches, SIAM/ASA Journal on Uncertainty Quantification 5 (1), 334-352, 2017 3. Litvinenko A., Matthies H.G., El-Moselhy T.A. (2013) Sampling and Low-Rank Tensor Approximation of the Response Surface. In: Dick J., Kuo F., Peters G., Sloan I. (eds) Monte Carlo and Quasi-Monte Carlo Methods 2012. Springer Proceedings in Mathematics & Statistics, vol 65. Springer, Berlin, Heidelber 4. A. Litvinenko, Application of hierarchical matrices for solving multiscale problems, Dissertation, Leipzig University, Germany, http://publications.rwth-aachen.de/record/754296/files/754296.pdf 5. A Litvinenko, R Kriemann, MG Genton, Y Sun, DE Keyes, HLIBCov: Parallel hierarchical matrix approximation of large covariance matrices and likelihoods with applications in parameter identification MethodsX 7, 100600, 2020 6. A Litvinenko, D Logashenko, R Tempone, G Wittum, D Keyes, Solution of the 3D density-driven groundwater flow problem with uncertain porosity and permeability. Int J Geomath 11, 10 (2020). https://doi.org/10.1007/s13137-020-0147-1 7. A Litvinenko, Y Sun, MG Genton, DE Keyes, Likelihood approximation with hierarchical matrices for large spatial datasets, Computational Statistics & Data Analysis 137, 115-132, 2019 8. S. Dolgov, A. Litvinenko, D.Liu, KRIGING IN TENSOR TRAIN DATA FORMAT Conf. Proceedings, 3rd International Conference on Uncertainty Quantification in Computational Sciences and Engineering, https://files.eccomasproceedia.org/papers/e-books/uncecomp_2019.pdf pp 309-329, 2019 9. A. Litvinenko, D. Keyes, V. Khoromskaia, B.N. Khoromskij, H. G. Matthies, Tucker tensor analysis of Mat´ern functions in spatial statistics, J. Computational Methods in Applied Mathematics, Vol. 19, Issue 1, pp 101-122, 2019, De Gruyter 10. HG Matthies, E Zander, BV Rosic, A Litvinenko, Parameter estimation via conditional expectation: a Bayesian inversion, Advanced modeling and simulation in engineering sciences 3 (1), 1-21, 2016 11. M Espig, W Hackbusch, A Litvinenko, HG Matthies, E Zander, Post-Processing of High-Dimensional Data, arXiv:1906.05669, 2019 33/33
- 35.