Algorithms for high performance computing applied to simulation of ensembles of nanoparticles

Abstract

(English) The objective of this thesis is to increase the computational efficiency of integral equation based methods for solving the electromagnetic scattering problem by incorporating random computations. Randomization allows to ease the computational burden by replacing costly exact operations by inexpensive approximate ones. Besides, it minimizes the need for inter-processor communications, leading to highly scalable algorithms. We select electromagnetic scattering in clusters of plasmonic nanoparticles as our case of study. Randomized methods are suitable for these kind of problems since their collective behavior can often be efficiently characterized by statistical means. The electromagnetic scattering problem in the integral form is discretized by the Method of Moments, which consists on expanding the unkown of the problem as a linear combination of basis functions and projecting both members of the equation onto a set of weighting functions. The result is a linear system of algebraic equations that approximates the original functional equation. This linear system is usually solved by iterative methods. Also, computing and manipulating the system matrix is unfeasible for most realistic problems. This is resolved by compression methods that take advantage of the reduced number of degrees of freedom of the electromagnetic interaction between distant subdomains of the scattering object. In this work, we propose a method for solving the electromagnetic scattering problem based on the Born-Neumann series. This formalism expresses the inverse of the scattering operator as a power series of the operator itself. Due to the structure of this expansion, Monte Carlo methods can be employed to accelerate the computation of the solution of the scattering problem. This thesis presents contributions on the following two topics: First, the relation of the Born-Neumann series with semiiteratiive methods, Krylov subspace methods and conformal mapping. Both Krylov and semiiterative methods approximate the inverse of a linear operator as a linear combination of powers of the operator itself. The difference lies on the way of finding the coefficients of such linear combination: whereas Krylov subspace methods actively find these coefficients by minimizing a certain residual, semiiterative methods compute them by applying a fixed transformation based on assumptions about the operator spectrum. We analyze the relation between both families of methods and study their connection with conformal maps. Second, the convergence properties of the Ulam-Neumann Monte Carlo algorithm. This method approximates the solution of a linear equation by performing random sampling on the terms of the Born-Neumann expansion. We study the convergence properties of this method when submatrices are sampled instead of single matrix elements, and propose an alterantive sampling scheme based on independently approximating each term of the Born-Neumann series that converges on a wider range of situations. In the last part of the work, we implement and test a randomized compression method for the impedance matrix obtained by the Method of Moments. The inherent parallelism of this method allows to exploit the dense linear algebra capabilities of high throughput processors as GPUs.

(Català) L'objectiu d'aquesta tesi és incrementar l'eficiència computacional dels mètodes basats en equació integral per resoldre el problema de la dispersió electromagnètica incorporant-hi càlculs aleatoris. L'aleatorització permet alleugerir la càrrega computacional substituint costoses operacions exactes per d'altres aproximades amb cost inferior. A més, minimitza la necessitat de comunicacions entre processadors, la qual cosa resulta en algorismes altament escalables. Escollim el problema de la dispersió electromagnètica en agrupacions de nanopartícules plasmòniques com a cas d'estudi. Els mètodes estocàstics són adequats per a aquesta mena de problema ja que el seu comportamnet col·lectiu sovint pot ésser caracteritzat per mitjans estadístics. El problema de la dispersió electromagnètica en forma integral es discretitza amb el Mètode dels Moments, que consisteix a expandir la incògnita del problema com a una combinació lineal de funcions base i projectar ambdós membres de l'equació sobre un conjunt de funcions pes. El resultat és un sistema lineal d'equacions algebraiques que aproxima l'equació funcional original. Aquest sistema lineal es resol habitualment mitjançant mètodes iteratius. D'altra banda, calcular i manipular la matriu del sistema és inviable en la majoria de casos realistes. Aquest problema se soluciona amb mètodes de compressió que aprofiten el reduït nombre de graus de llibertat de les interaccions electromagnètiques entre subdominis llunyans de l'objecte dispersiu. En aquest treball proposem un mètode per a resoldre el problema de la dispersió electromagnètica basat en la sèrie de Born-Neumann. Aquest formalisme expressa la inversa de l'operador de dispersió com a una sèrie de potències del propi operador. Degut a l'estructura d'aquesta expansió, els mètodes de Montecarlo poden ser emprats per a accelerar el càlcul de la solució del problema de la dispersió. Aquesta tesi presenta contribucions en els dos àmbits següents: En primer lloc, la relació entre la sèrie de Bonr-Neumann i els mètodes semiiteratius, els mètodes de subespai de Krylov i les transformacions conformes. Tant els mètodes de Krylov com els mètodes semiiterarius aproximen la inversa d'un operador lineal com a combinació lineal de potències del propi operador. La diferència rau en la forma de trobar els coeficients d'aquesta combinació lineal: mentre que els mètodes de subespai de Krylov busquen activament aquest coeficients minimitzant un cert residu, els mètodes semiiteratius els calculen aplicant una transformació fixa basada en assumpcions sobre l'espectre de l'operador. Analitzem la relació entre ambdues famílies de mètodes i estudiem les seves connexions amb les transformacions conformes. En segon lloc, les propietats de convergència del mètode de Montecarlo d'Ulam-Neumann. Aquest mètode aproxima la solució d'un sistema d'equacions lineals executant un mostreig aleatori sobre els termes de l'expansió de Born-Neumann. Estudiem les propietats de convergència d'aquest mètode quan es mostregen submatrius en lloc d'elements individuals de la matriu, i proposem un esquema de mostreig alternatiu basat en aproximar independentment cada terme de la sèrie de Born-Neumann que convergeix en un ventall més ample de situacions. A la darrera part del treball, implementem i provem un mètode aleatori de compressió de la matriu d'impedàncies obtinguda pel Mètode dels Moments. El paral·lelisme inherent d'aquest mètode permet explotar les capacitats per a l'àlgebra lineal de processadors d'alt rendiment com les GPUs.

(Español) El objetivo de esta tesis es incrementar la eficiencia computacional de los mètodos de ecuación integral para resolver el problema de la dispersión electromagnètica incorporando cálculos aleatorios. La aleatorización permite aligerar la carga computacional sustituyendo costosas operacions exactas por otras aproximadas de menor coste. Además, minimiza la necesidad de comunicaciones entre procesadores, dando lugar a algoritmos altamente escalables. Elegimos el problema de la dispersión electromagnètica en agrupaciones de nanoportículas plasmónicas como caso de estutdio. Los métodos estocásticos son adecuados para este tipo de problemas porque a menudo pueden ser caracterizados por medios estadísticos. El problema de la dispersión electromagnética en su forma integral se discretiza con el Método de los Momentos, que consiste en expandir la incógnita del problema como una combinación de funciones base y proyectar ambos lados de la ecuación sobre un conjunto de funciones peso. El resultado es un sistema lineal de ecuaciones algebraicas que aproxima la ecuación funcional original. Este sistema lineal se resuelve habitualmente con métodos iterativos. Por otro lado, calcular y manipular la matriz del sistema is inviable para la mayoría de problemas realistas. Esto se soluciona mediante el uso de métodos de compresión que aprovechan el reducido número de grados de libertad de la interacción electromagnética entre subdominios distantes del objeto. En este trabajo, proponemos un método para resolver el problema de la dispersión electromagnética basado en la serie de Born-Neumann. Este formalismo expresa la inversa del operador de dispersión como una serie de potencias del propio operador. Debido a la estructura de esta expansión, los métodos de Montecarlo pueden emplearse para acelerar el cálculo de la solución del problema de la dispersión. Esta tesis presenta contribuciones en los siguientes dos ámbitos: En primer lugar, la relación entre la serie de Born-Neumann, los métodos semiiterativos, los métodos de subespacio de Krylov y las transformaciones conformes. Tanto los métodos de Krylov como los semiiterativos aproximan la inversa de un operador lineal como una combinación lineal de potencias del propio operador. La diferencia estriba en la forma de encontrar los coeficientes de tal combinación lineal: mientras que los métodos de subespacio de Krylov buscan activamente estos coeficientes minimizando un cierto residuo, los métodos semiiterativos los obtienen aplicando una transformación fija basada en asunciones acerca del espectro del operador. Analizamos la relación entre ambas familas de métodos y estudiamos su conexión con las transformaciones conformes. En segundo lugar, las propiedades de convergencia del método de Montecarlo de Ulam-Neumann. Este método aproxima la solución de la ecuación mediente un muestreo aleatorio sobre los términos de la expansión de Born-Neuman. Estudiamos además las propiedades de convergencia de este método cuando se muestrean submatrices en lugar de elementos individuales de la matriz, y proponemos un esquema alternativo de muestreo basado en aproximar independientemente cada término de la serie de Born-Neumann que converge en un abanico más amplio de situaciones. En la última parte de este trabajo, implementamos y probamos un método de compresión estocástico para la matriz de impedancias obtenida con el Método de los Momentos. El paralelismo inherente de este método permite explotar las capacidades de procesadores de alto rendimiento como las GPUs.