Alexis Gobé

Professeur Christophe Geuzaine

Méthodes Galerkin discontinues pour la simulation de problèmes multiéchelles en nanophotonique et applications au piégeage de la lumière dans des cellules solaires

Keywords: Équations de Maxwell, Méthode de Galerkine discontinue, Matériaux dispersifs, Photovoltaique, Piégeage de la lumière, Nanophotonique, Maillage hybride, Maillage non conforme, Méthode d Optimisation, Calcul haute performance Maxwell's equations, Discontinuous Galerkin method, Dispersive materials, Photovoltaics, Light trapping, Nanophotonics, Hybrid mesh, Non conform mesh, High order method, Multiscale method, optimization, High performance computing

Titre: Méthodes Galerkin discontinues pour la simulation de problèmes multiéchelles en nanophotonique et applications au piégeage de la lumière dans des cellules solaires Résumé: L'objectif de cette thèse est l'étude numérique du piégeage de la lumière dans des cellules solaires nanostructurées.

Le changement climatique est devenu une problématique majeure nécessitant une transition énergétique à court terme. Dans ce contexte, l'énergie solaire semble être une source énergétique idéale. Cette ressource est à la fois scalable à l'échelle planétaire et écologique. Afin de maximiser sa pénétration, des travaux visant à augmenter la quantité de lumière absorbée et à réduire les coûts liés à la conception des cellules sont nécessaires. Le piégeage de la lumière est une stratégie qui permet d'atteindre ces deux objectifs. Son principe consiste à utiliser des texturations nanométriques afin de focaliser la lumière dans les couches de semi-conducteur absorbantes. Dans ce travail, la méthode de Galerkine Discontinue en Domaine Temporel (DGTD) est introduite. Deux développements méthodologiques majeurs, permettant de mieux prendre en compte les caractéristiques des cellules solaires, sont présentés. Tout d'abord, l'utilisation d'un ordre d'approximation local est proposé , basé sur une stratégie de répartition particulière de l'ordre. Le deuxième développement est l'utilisation de maillage hybride mixant ses élements hexahédriques structurés et tétrahédriques non structurés. Des cas réalistes de cellules solaires issus de la littérature et de collaborations avec des physiciens du domaine du photovoltaïque permettent d'illustrer l' apport de ces développements. Un cas d'optimisation inverse de réseau de diffraction dans une cellule solaire est également présenté en couplant le solveur numérique avec un algorithme d'optimisation bayésienne. De plus, une étude approfondie des performances du solveur a également été réalisée avec des modifications méthodologiques pour contrer les problèmes de répartition de charge. Enfin, une méthode plus prospective, la méthode Multiéchelle Hybride-Mixte (MHM) spécialisée dans la résolution de problème très fortement multiéchelle est introduite. Un schéma en temps multiéchelle est présenté et sa stabilité prouvée.

Si une partie de ma thèse devait être écrite en français, c'est bien celle là. Au moment où j'écris ces mots, je suis exactement à un mois du point finale de cette aventure. Quelle aventure, quelle difficile et éprouvante aventure, quelle agréable aventure. Je remercie tout d'abord mon boss et directeur de thèse Stéphane, qui a toujours été là pour moi. On m'a souvent parlé de relations thésard/directeur de thèse difficiles, ou quasiment inexistantes. Pour moi ça a été tout l'inverse, et je ne changerai absolument rien à ces 3.7 années à avoir travaillé avec toi. Claire, merci pour tout. Tu as été indispensable durant ces 3 années. Tu as notamment réussi à me faire comprendre que la recherche, ce n'était pas seulement trouver, et que trouver quelque chose qui ne marche pas trop, c'était avoir trouver quelque chose! Merci à Stéphanie et Christophe pour avoir accepter de rapporter sur ce manuscrit, vos remarques m'ont été particulièrement utile pour sa finalisation. Merci à Amélie, Boniface, Claire et Stéphane d'avoir accepter d'être éxaminateurs. J'ai particulièrement apprécié les discussions scientifiques que nous avons eues avant, pendant et après la soutenance. Jonathan, merci pour tout. Cette thèse n'aurait clairement pas la même saveur sans toi et ton goût pour l'esthétisme et pour tikz. Épargne moi tes "voleur" par ci, "voleur" par là en voyant mes planches de soutenance. . . ! Niko, danke! Sehr gut café! C'était super d'être doctorant à tes côtés, toujours de supers discussions, que ce soit des discussions scientifiques ou des discussions sur les chars de combat. . . Docteur Léger, mon voisin mouginois du haut, merci beaucoup pour ce passage de témoin avec la MHM lors de tes derniers mois dans l'équipe. Je te serai éternellement reconnaissant pour ton accueil à Petropolis. C'est sûr que mon séjour au LNCC aurait été très très très différent sans toi. Ludo merci pour ton aide pendant TECSER, c'est grâce à toi que mon stage c'est si bien passé, ce qui m'a permis de commencer cette thèse. Merci Alexandra, surtout pour cette séance de sport ensemble (j'ai encore mal au jambe). Merci George, que ce soit pour le taxi, pour les burgers, pour les pâtisseries libanaises. . . Le prochain c'est toi! Je vais essayer de remercier également tous les Nachos/Atlantis que j'ai vu passer et avec qui j'ai passé de bons moments (Mahmoud, Nora, Vincent, Mohamed, Théo, Kun, Plum, Omid, Hao, Guillaume, Patrick, Théo, Yves, Stéphane D, Montserrat, Zakaria, Massimiliano, Thibault, Cédric, Leandro, Vincent, . . . ), mes collègues d'autres équipes (Mickael, Pascal, Radu, Quentin, Milica, Nicholas, Lorenzo, Émilie,...), et mes rencontres d'autre pays (Wesley, Frederic, Tadeu, ...). Je remercie mes parents, sans qui je n'aurai été capable de rien. Merci d'avoir toujours été là, tout le temps depuis maintenant 28 ans. Je n'aurai jamais la possibilité de vous exprimer toute la reconnaissance que j'ai pour vous. Mémé, merci pour tout. J'ai tellement de chance de vous avoir. Quand j'y pense c'est difficile d'être à 1000 kilomètres de vous. Je vous dédis ces travaux. Merci également à ma famille au sens plus large, oncles, tantes, Audrey, Valérie, Hervé, . . . Merci à mes deux grands amis, Clément et Maxence, pour être là depuis bientôt 20 ans. Je vous remercie tout particulièrement pour tous ces soirs où j'étais seul la semaine, et où nos soirées sur le pc à jouer et à parler ont permis à ces 3 années de passer bien plus vite que prévu. Merci également aux potes de la fac (Amélie, Flore, Léa, Pauline, ...) et à tous les autres que j'oublie. Enfin, j'ai gardé le meilleur pour la fin, enfin la meilleure. Docteur Dupré. . . enfin pas encore mais c'est une question de jour maintenant! Sans toi, rien de tout ça ne se serait passé. Sans toi, c'est certain que je n'aurais jamais fait de thèse, ni de master de modélisation mathématique, ni de master de mathématiques d'ailleurs. C'est vraiment une chance de t'avoir eu sur ma route ce mois de mars 2012. On dira qu'on a survécu à ces 3 années à devoir attendre le vendredi soir pour se voir. Que c'était difficile. Merci d'avoir fait toutes ces concessions pour moi, pour nous... et surtout d'avoir fait exploser ton indice carbone, c'est plutôt cocasse quand on a fait une thèse sur une energie renouvelable.

Introduction

Mathematical models can describe a wide variety of phenomena, which can be physical, biological, economic, demographic, geological for example. There are a considerable number of models, from different families, sometimes established well before the birth of the first computers. Some are stochastic, e.g. with consideration of randomness, while others are deterministic, e.g. in the form of an ordinary differential equation (ODE) or a system of partial differential equations (PDE). Those models can be used to simulate realistic phenomena, which can thus play a predictive or explicative role. The study and exploitation of these models have always been essential for engineers. The importance of numerical modeling has become a more significant challenge with the increase in computer system capabilities, and in particular, with the advent of parallel computers. This was particularly the case for electromagnetism. The physicist James Clerk Maxwell laid the foundations of modern electromagnetism in 1865 in [START_REF] Maxwell | A Dynamical Theory of the Electromagnetic Field Email alerting service[END_REF]. For the first time, he formulated the classical theory of electromagneticm by bringing together electricity, magnetism, and light as various manifestations of the same phenomenon. Maxwell's equations have been studied for many decades, for different purposes. As for most of the computational models, computational electromagnetism can be used for validation and prototyping in a goal-oriented engineering setup. We can find electromagnetism applications everywhere and every day. Many of those are part of our daily life, such as wireless communications, electromagnetic braking, optical fibers, medical imaging, magnetic lift train, induction charging of batteries, or even in our computer's hard drives.

Photonics is one of the fields modeled by Maxwell's equations. It is the physical science of light, which is an electromagnetic radiation. One of the subfields of photonics, which is at the heart of this thesis work, is called nanophotonics.

Nanophotonics and photovoltaics Nanophotonics is the study of light and its interactions with matter at nanometric scales (or nanoscale). Nanophotonics incurs phenomena occurring at distances smaller or equivalent to the wavelength of the light. In the visible range of the electromagnetic spectrum (see figure 1.1), the wavelength is between 400 nm (which is the ultraviolet limit) and 700 nm (which is the infrared limit). When the studied structure is smaller or approaching this particular scale, the light-matter interactions are unveiling new phenomena, and general optical physics is no longer appropriate. At this range, it is important to take into account the wavelength dependency of the speed of light. This phenomenon is called dispersion, and several models aim to simulate this material specificity, as the Drude model [START_REF] Drude | Zur elektronentheorie der metalle[END_REF], and the Generalized dispersive model [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the fdtd method[END_REF][START_REF] Viquerat | Fitting experimental dispersion data with a simulated annealing method for nano-optics applications[END_REF]. Over the past decades, nanophotonics research has led to astonishing technologies. To name some of the most important, the metasurfaces, which are thin sheet metal surfaces that allow the reflected and transmitted waves to be modulated in phase [GCA + 17]. Thanks to optothermal heating [SSD + 14], nanoparticles are now used in medicine in order to kill cancer cells. Advanced metasurfaces based on nanocube are now being used for nanoscale holography [START_REF] Genevet | Holographic optical metasurfaces: a review of current progress[END_REF] and very precise molecule sensing [PCT + 16]. Researchers have investigated a variety of nanophotonic methods in order to concentrate light in the optimal areas within a solar cell [START_REF] Ferry | Design Considerations for Plasmonic Photovoltaics[END_REF]. Improvements in fabrication processes have permitted to reach high performances up to 44.7% [DGB + 14].

In the context of climate change and energy transition, the development of renewable energies is essential. In the world, they represent a little more than 20% of the energy mix. If we put aside hydropower (which is a particular renewable energy because there is no intermittency), the first renewable resource in the world is wind power (9% in 2017) followed by solar energy (4% in 2017). However, according to the International Energy Agency (IEA), renewable energies will grow most rapidly in the electricity sector, providing nearly 30% of electricity demand by 2023, as shown. During this period, renewable energies are expected to account for more than 70% of global electricity generation growth, led by solar photovoltaics, followed by wind, hydro, and bioenergy. In particular, solar photovoltaics is expected to significantly increase all over the world, as shown in figure 1.2. Also, the IEA1 predicts that photovoltaics could produce up to 16% of the world's electricity by 2050, thus becoming the first renewable resource.

With this in mind, a great amount of research is being carried out to improve efficiency and reduce the fabrication cost of solar cells. "To support global research efforts to produce photovoltaic modules with a >30% energy conversion efficiency for a <30 c$/Wp price by 2030: [30/30/30]" is the roadmap put forth by representatives of the leading solar photovoltaics research institute in France (IPVF2 ). Achieving this objective requires several major innovations. Already promising technologies, such as multijunction thin-films and crystalline silicon solar cells, are good candidates, but for the moment, they suffer from prohibitive manufacturing cost. The development of new experimental concepts based on more complicated physics, such as hot carriers are also considered. Finally, the use of nanophotonics in order to trap light inside the solar cell is an appealing possibility, mostly driven by the decrease in the thickness of solar cells to reduce materials costs [START_REF] Mokkapati | Nanophotonic light trapping in solar cells[END_REF]. Nanoscale light trapping allows the absorption to exceed the limits attainable with geometrical optics.

Light trapping Light trapping consists in capturing photons for use in several applications, such as solar cells. Light enters the different structures by refraction or scattering. Some light trapping structures allow incident light to be guided into the absorbent layer of a solar cell. The purpose is to keep the light trapped inside the device so that the absorption capacity can be improved beyond the absorption of a single pass of light through the device (see figure 1.3).

Different types of nanophotonic structures have been reported in the literature, such as gratings [WYL + 12], nanowires [START_REF] Eyderman | Near perfect solar absorption in ultra-thin-film {GaAs} phtonic crystals[END_REF] and random scattering [START_REF] Lockau | Optical modeling of thin film silicon solar cells with random and periodic light management textures[END_REF]. These structures are very difficult to choose a priori (before any measures or simulations). Physicists generally need to go through optimization phase, both parametric and geometric, to find the design that will lead to the best efficiency. Computational electromagnetism Nowadays, various numerical methods are well known and used to solve Maxwell's PDE system. These methods are especially useful and needed when the application to simulate cannot be handled with analytical solutions. Numerical methods make it possible to solve complex problems, quickly validate a corrective action, solve complex problems by splitting them into subproblems, optimize the sizing of installations, access non-measurable quantities, better understand the physical phenomena involved, but above all, reduce costs by solving problems before they occur. However, most of the realistic electromagnetic problems cannot be solved analytically. The complexity of the underlying PDE systems does not allow for analytical solutions except in a few simple geometrical settings. Therefore, the use of numerical methods is the only way to study complex and real-life problems.

The full-wave Maxwell's equations can be formulated in two settings: the time-domain (TD) and the frequency-domain (FD) formulations. Several numerical methods have been devised to solve both formulations. As often in numerical modeling, the choice of a method is driven by the application.

Frequency-domain methods Frequency-domain methods, also known as time-harmonic methods, are solving Maxwell's equations in frequency-domain. This means that the time variation is supposed to be harmonic. Frequency-domain methods are used when the considered propagation problem is monochromatic. The Finite Element Frequency-Domain (FEFD) is used to solve the frequency-domain boundary valued electromagnetic problems by using a variational form [Can]. It can be used with two-and three-dimensional canonical elements of differing shapes, allowing for a highly accurate discretization of the solution domain. The FEFD is often used in the frequency domain for computing the frequency field distribution in complex, closed regions such as cavities and waveguides. However, the solution domain must be truncated, making the FEFD unsuitable for radiation or scattering problems. The Method Of Moments (MOM) is a technique used to solve electromagnetic boundary or volume integral equations in the frequency-domain [BVF + 17]. Because the electromagnetic sources are the quantities of interest, the MOM is very useful in solving radiation and scattering problems. Rigorous Coupled-Wave Analysis (RCWA) method is a semianalytical method. It is a Fourier-space method, so devices and fields are represented as a sum of spatial harmonics [START_REF] Moharam | Rigorous coupled-wave analysis of planar-grating diffraction[END_REF]. This method is used to improve diffraction structures for high efficiency solar cells. For not so complex 3D structure, RCWA is one of the reference methods. However, when the structure studied is varying in the propagation direction, RCWA needs to discretize it in a set of uniform (in the propagation direction) layer.

Time-domain methods Unlike the frequency-domain methods, time-domain methods with explicit time integration may not require to solve a large linear system of equations. Besides, it is possible to compute the solution over a broadband of wavelengths in a single simulation. For many years, the reference method for solving nanophotonic and PV broadband problems has been the Finite Difference Time-domain (FDTD) method or Yee's method [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF]- [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]- [START_REF] Gallinet | Numerical methods for nanophotonics: Standard problems and future challenges[END_REF]. The main reasons for the popularity of this method are its ease of implementation and its computational efficiency. However, to achieve the latter properties, the FDTD relies on a cartesian structured discretization of the computational domain [START_REF] Holland | Finite-Difference Analysis of EMP Coupling to Thin Struts and Wires[END_REF]. When the underlying physics involves local rounded details, the structured grid becomes a significant limitation, which is manageable but makes the method more complex. Being naturally constructed for unstructured grids, the Finite Element method (FEM) is well suited for numerical solution of PDEs in complex geometries [START_REF] Cangellaris | Point-Matched Time Domain Finite Element Methods for Electromagnetic Radiation and Scattering[END_REF]. Although initially developed for elliptic equations (boundary value problems), FEM have been further extended to hyperbolic equations (initial boundary value problems). In the time-domain and for the classic FEM, one generally needs to solve the sparse linear system associated to the mass matrix at each time-step, which makes the method very expensive and rarely used in practice. The Finite Volume method (FVM) relies on the evaluation of integrals over volume elements and facets of the mesh. This allows great geometrical flexibility, which is required for most electromagnetic problems [START_REF] Fumeaux | Developments of finitevolume techniques for electromagnetic modeling in unstructured meshes[END_REF]. However, this method lends itself to a piecewise constant approximation of the field and is thus facing accuracy issues. The Boundary Element Method (BEM) is an important alternative to the prevailing domain methods, such as the FDTD or the FETD. The basis of the method is to transform the original system of PDEs into an equivalent integral equation (or system) using the corresponding Green's representation formula (direct method), or in terms of a continuous distribution of singular solutions of the PDE over the boundaries of the problem (indirect method) [START_REF] Huybrechs | An efficient implementation of boundary element methods for computationally expensive Green's functions[END_REF]. Finally, about twenty years ago, the Discontinuous Galerkin for time-domain Maxwell's equations (DGTD) was introduced in [START_REF] Remaki | Une méthode de Galerkin discontinu pour la résolution des équations de Maxwell en milieu héterogene[END_REF].

Discontinuous Galerkin Time-Domain Among the previously presented methods, the Discontinuous Galerkin Time-Domain (DGTD) method has gradually emerged as a interesting alternative to the FDTD method. It is a discontinuous finite element type method that relies on a high-order interpolation of the electromagnetic fields within each cell of an unstructured mesh. Time integration can be achieved using an explicit scheme and, as a result of the discontinuity of the approximation, no global mass matrix inversion is required to advance the solution at each time-step, as it is the case with FETD methods working on unstructured meshes. However, the use of an explicit scheme restricts the time-step used, respectively to the Courant-Friedrichs-Lewy (CFL) condition. The resulting time-domain solver is well adapted to massively parallel computing. DGTD can also easily handle complex meshes thanks to the natural discontinuity of the solution at the cell's interface. The grid may be a conforming unstructured finite element mesh, a non-conforming one, or even a hybrid mesh made of various elements (tetrahedra, hexahedra) without the need of transition elements (prisms, pyramids) in contrast with continuous finite element methods. The DGTD method has been proven to work well with highly locally refined meshes [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF]. This property makes this method more suitable for designing a hp-adaptive solution strategy (i.e. where the characteristic mesh size h and the interpolation degree p changes locally wherever it is needed). It is also flexible concerning the time integration scheme. The discontinuous Galerkin spatial discretization can be combined with any explicit or even globally or locally implicit time-stepping scheme.

Multiscale problems Many applications considered in nanophotonics and photovoltaics involve highly heterogeneous media or/and nanotextures embedded in larger domains. Such physical problems are particularly challenging for numerical methods. Indeed, whenever multiple spatial or temporal scales play an essential role in the problem of interest, several numerical difficulties must be addressed. In particular, discretization methods have to use very fine meshes to represent the smallest geometrical details. More specifically, the main difficulty from a numerical modeling point of view is to find parameters for constitutive relationships from material models that are valid at a smaller scale. Indeed, few standard algorithms are efficient, accurate, and robust enough to solve Maxwell's macroscopic equations when calculations involve multiple scales. In recent years, discretization methods suitable for multiscale problems have been studied extensively, such as Discontinuous Galerkin (DG) with strong heterogeneity of the electric permittivity as well as very small and curved geometries [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF][START_REF] Chen | Discontinuous Galerkin time-domain methods for multiscale electromagnetic simulations: A review[END_REF]. While these two methods are quite general, homogenization techniques such as Heterogeneous Multiscale (HMM) method [START_REF] Henning | A new heterogeneous multiscale method for time-harmonic Maxwell's equations[END_REF], or the use of multiscale basis functions as in the Hybrid High-Order (HHO) method [START_REF] Di Pietro | An introduction to Hybrid High-Order methods[END_REF] help dealing with both scales. The recently introduced MHM method, originally presented in [START_REF] Araya | Multiscale Hybrid-Mixed Method[END_REF] for elasticity and recently extended for Maxwell equations in time-domain [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF], seems to be a promising approach, both in terms of flexibility and performance thanks to its divide-and-conquer framework totally suitable for massively parallel computer architectures.

In this thesis, we want to show the interest of the DGTD method for the simulation of photovoltaic problems. In particular, the trapping of light in solar cells involves both geometric and physical problems, which allows the use of possibilities of the DGTD method, which is rarely used. The MHM method, still very recent for Maxwell's equation in the time domain, seems very promising for solving photovoltaic problems on a larger scale.

Contributions and thesis outline

The manuscript is structured as follow:

• In chapter 2, we introduce the photovoltaic topic, which is nowadays an essential concern in this context of energy transition. We present the main concepts of PV cells, and especially the optical part, which explains how photons are absorbed in the solar cell. As optics is directly related to electromagnetism, we recall the system of Maxwell's equations, as well as its extension to deal with the peculiarities of nanophotonics. Then, we present the spatial discretization of Maxwell's equations using a DG method. We briefly present the time-stepping strategy chosen for the time discretization. Additional details about domain truncation and source definition are also presented.

• In chapter 3, we present extensions to the standard Discontinuous Galerkin Time-domain (DGTD) method that we have studied in the context of the simulation of complex multiscale PV structures. We propose two numerical strategies, which are well-known for being especially well suited to DG methods, in order to deal with the specificities of solar cell structures, in particular, the periodic nanostructured stacking layer. A p-adaptivity strategy is proposed and compared on a realistic PV cell with the one initially present in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF]. We show nice improvements, both in terms of accuracy and computational time. We also present the DGTD method on hybrid structured-unstructured meshes, which allows us to improve the computational efficiency of realistic cases further.

• Chapter 7 is dedicated to the MHM for time-domain Maxwell's equations. The MHM is a twolevel method, which requires the use of a separate numerical method at the second level. For the latter, we choose to use the DGTD method. This MHM method was initially introduced in [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF], and is extended here with a multiscale time-stepping scheme in order to get the most out of the MHM. The stability is theoretically proved, and 2D results are presented. For the first time, we show some numerical validations (energy preservation, convergence study), and then we show the capability of the method on a realistic nanophotonic case.

• Chapter 4 is concerned with realistic PV applications. We start with an introduction of stateof-the-art methods to deal with nanophotonic cases. A brief explanation of each method's advantages and drawbacks is proposed, and DG is placed in relation to them. Then, we define the main figures of merit used in photovoltaics and how they are computed in the DGTD framework. Finally, we present how realistic measures on solar cell nanotextures can be handled to generate high fidelity meshes. Two strategies are proposed in order to obtain biperiodic structures of complex shapes.

• In chapter 5, we concentrate on the optimization of PV cells. We investigate the coupling of an optimization method called efficient global optimization (EGO), with our DGTD solver. This optimization method relies on a metamodelisation. The metamodels are used to mimic the figure of merit, at a very cheap cost compared to the full-wave DGTD solver. Two realistic PV optimization problems are studied using the EGO methods. Efficiency and performance are presented, in particular with the parallel EGO (q-EGO), which allows parallel evaluation of the DGTD solver.

• Finally, in chapter 6, we study the high-performance computing aspects of our DGTD solver. We investigate the effect of traditional finite element methods parallelization, which is based on a mesh partition. We show that we have to face load balancing issues, mostly due to local computations inherent to realistic PV simulations. Advanced strategies are presented in order to handle really complex configurations and evaluated on realistic PV problems.

In the context of climate change, the global energy market is turning increasingly into green energies since the end of last century. Renewable energies are scalable to world energy demand. Wind, solar, or hydro energies are the leading candidates to replace fossil energies like nuclear, coal, and gas. Several studies have aimed at estimating the potential of each renewable energy technology, and solar energy has been identified as the favorite [START_REF]Global trends in renewable energy investment[END_REF]. Solar energy is an abundant, clean, and available source of power. In 2017, more than 37% of power generating capacity installed was solar installations (see figure 2.1). This is nearly twice as much capacity as wind power [START_REF]Renewables[END_REF]. With a potential of 1575 exajoule (EJ) per year, while world energy consumption was 606 EJ in 2015 [START_REF] Goldemberg | World energy assessment. energy and the challenge of sustainability[END_REF], solar energy could ensure the world's energy demand in the near future. The physics of solar cells is based on the photovoltaic effect, discovered by Edmond Becquerel [START_REF] Becquerel | Memoire sur les effets electriques produits sous l'influence des rayons solaires[END_REF] in the 19th century. The main idea behind a solar cell is to convert light into electrical power. Historically, space programs were the largest supporters of PV technology since photovoltaics was the best energy source for their satellites. The industry has since grown, and one can have seen PV systems almost everywhere in our daily life. For example, it is used to power cars, houses and commercial buildings.

PV cells are embedded in solar panels, which are also part of the so-called PV system (see figure 2.2). The generation of solar energy requires optimizing each part of a photovoltaic system, In order to continue this massive spread of PV around the world, it is necessary to keep improving the efficiency of solar cells while lowering manufacturing and material prices. One particularly well suited approach is to use nanophotonics and light trapping in order to maximize the absorption of the light inside the solar cell. This strategy allows to reduce the thickness of solar cells, which naturally decrease materials costs [START_REF] Mokkapati | Nanophotonic light trapping in solar cells[END_REF]. Most of the simulation of nanophotonics phenomena relies on the resolution of the system of Maxwell's equations.

Apart from basic devices (like rods or spheres) where analytical solutions are available (mostly thanks to Mie theories), realistic configurations are most of the time impossible to solve by hand. For that reason, numerical methods which solve Maxwell's equations have received considerable interest in the last decades. The family of DGTD methods seems promising. DGTD methods rely on a weak formulation and work very well on unstructured meshes. Their main specificity comes from the discontinuous approximation space used that decouples the neighbored mesh elements. Such an approach leads to discontinuities at the cell interfaces and hence to a global discontinuous solution. Moreover, the decoupling of the mesh elements makes the mass matrix block-diagonal, which can therefore be easily invert. This makes DGTD well-suited for explicit time integration schemes.

The first part of this chapter (section 2.1) is an introduction to photovoltaics. Since this thesis is not specifically about the physical aspects of solar light energy harvesting, we intentionally keep the discussion at a general level. This section includes a brief overview of the most used solar cell technologies. Although the numerical modeling of a full solar cell device required to deal with multiple physical models, i.e., it is a coupled multiphysics problem, our main concern in this thesis is the characterization of solar light interaction with a solar cell structures. The optical part of photovoltaics is modeled thanks to the system of Maxwell equations. Theses equations are recalled in section 2.2. Since we adopt a time-domain setting, we discuss in more detail the modeling of material physical dispersion and its role in the present context. The third section of this chapter is dedicated to the numerical modeling of light with dispersive media and detail the high-order DGTD method that constitutes the starting point of this thesis work (section 2.3).

Photovoltaics

Single junction and multiple junction

Multi-junction solar cells are p-n junction solar cells made of different semiconductor materials. The p-stands for a deficiency of electrons and an abundance of electron holes, while n-stands for the abundance of electrons and few electron holes The p-n junction of each material produces an electrical current in response to different wavelengths of light. Thus, the use of various semiconductor materials allows the absorption of a larger range of wavelengths. The efficiency of converting sunlight from the cell into electrical energy is therefore improved. Traditional single-junction cells have a maximum theoretical efficiency of 33.7%, while an infinite number of junctions would have a limit efficiency of 86.8% (theoretically) under highly concentrated sunlight [START_REF] Green | Third generation photovoltaics: Ultra-high conversion efficiency at low cost[END_REF].

Semiconductor materials

Most of the PV cells are made of a stacking of layers of semiconductor materials. Semiconductor materials come from different parts of the periodic table, but they have some similarities. Their properties are related to their atomic characteristics and vary from one group to another. The atoms of a semiconductor can be materials from group IV of the periodic table, or a combination of material from group III and group V (called III-V semiconductors), or also combinations of material from group II and group VI (called II-VI semiconductors). Silicon is a group IV semiconductor material, and it is the most commonly used. One finds it in the basis of integrated circuit (IC) chips, which also make it the most advanced technology. Thus, most solar cells are made of silicon.

To (highly) schematize the PV phenomenon, the incoming light (photons) strikes the solar cell surface, and some of it is absorbed through the semiconductor, turning photons into electrons. The most famous PV cell is the single crystalline silicon, which consists of two differently doped silicon layers (p-n junction) [START_REF] Chapin | A new silicon p-n junction photocell for converting solar radiation into electrical power[END_REF] as represented in figure 2.3. A solar cell can be seen as a battery, pushing electrons from the p-to the n-layer.

Direct/indirect bandgap materials

The bandgap is the minimum amount of energy required for an electron to release itself from its bound state in order to participate in the conduction. The bandgap determines the amount of solar energy required for conduction, as well as the amount of energy produced. The bandgap of a semiconductor can be either a direct bandgap or an indirect bandgap. The difference between the two is based on the concepts of valence band and conduction band, as shown in figure 2. 4. In non-metals, the valence band is the highest range of electron energies in which they are normally present at absolute zero, while the conduction band is the lowest range of free electron states. In the case where the energy of the photon is higher than the bandgap, an electron-hole pair can be generated by exciting an electron of the valence in the conduction band. Any extra energy from the photon above the bandgap is lost due to thermalization.

A direct bandgap semiconductor is a semiconductor for which the maximum energy level of the valence band aligns with the minimum energy level of the conduction band, as shown in figure 2.4(a). Thus, a direct recombination takes place with the release of energy equal to the difference in energy between the recombinant particles. The probability of radiative recombination is high. The efficiency factor of such a semiconductor is higher than for the other bandgap. Thus, these semiconductors are always preferred over indirect bandgap semiconductors for the manufacture of optical sources. Among these semiconductors, one may find amorphous silicon and some III-V materials such as GaAs.

Within an indirect bandgap semiconductor, the maximum energy level of the valence band and the minimum energy level of the conduction band are misaligned, as shown in figure 2.4(b). In this case, due to a relative difference in the voltage pulse, the voltage pulse is first preserved by the release of energy, and only once the two momenta are aligned, a recombination occurs. The latter being accompanied by the release of energy. The probability of radiative recombination is relatively low, and the efficiency factor of such semiconductors is lower. These materials include crystalline silicon and Ge.

Efficiency and absorption

Absorption and efficiency of a PV system are two of the most essential properties of a solar cell. The efficiency of solar cells refers to the portion of sunlight energy which can be transformed by photovoltaics into electricity. This allows to determines the energy production of a PV system.

Absorption of light Three types of phenomena can occur when an incident photon impinges the surface of a semiconductor. The latter can be reflected from the top surface, it can be absorbed by the material, or it can be transmitted by the material. In the case of a solar cell, if the photon is reflected or transmitted, this will be considered a loss. Indeed, if the photon is not absorbed, it cannot produce energy. On the other hand, if it is absorbed, it can then excite an electron from the valence band to the conduction band. The response of the photon (absorption or transmission) is determined by its energy. If the photon does not have enough energy, it will not be able to excite the electron. The photon will therefore be transmitted. In order to be absorbed and excite the electron in the conduction band from the valence band, the incident photon must have sufficient energy. We can then define three groups according to the photon energy E ph and the energy of the semiconductor bandgap E bg :

• When the energy of the photons is lower than that of the material bandgap band E ph < E bg , the photons interact weakly with the semiconductor. In this case, they pass through it.

• When the energy of the photons is equal to that of the material bandgap E ph = E bg , the photons have just enough energy to create a pair of electron holes. In this case, photons are successfully absorbed.

• When the energy of the photons is higher than that of the material bandgap E ph > E bg , the photons are strongly absorbed. However, in the case of solar cells, the excess of energy is wasted due to the thermalization of electrons.

When a photon is absorbed, this results in the creation of a majority and a minority carrier. This production of charge carriers is the basis for PV energy production.

Absorption coefficient

The absorption coefficient can quantify the ability of a material to absorb photons. This coefficient determines how far the light can penetrate the material before it is fully absorbed. It is noted α and is defined as follows:

α = 4πk λ (2.1)
where k is the extinction coefficient and λ is the wavelength. When the absorption coefficient is low, this means that the material can weakly absorb the light. This coefficient therefore depends on the material, but also on the wavelength of the incident light. Figure 2.5 displays the absorption coefficient for three semiconductors: Gallium Arsenide (GaAs), Silicon (Si) and Amorphous Silicon (a-Si). The absorption coefficient is shown for the materials at 300°K depending on the wavelength of the light in vacuum. One can see a sharp edge in their absorption coefficient. Indeed, when the energy of the photons is lower than the bandgap as mentioned above, they cannot excite an electron from the valence band into the conduction band. Thus, electrons cannot be absorbed.

When the energy of the photons is close to the energy of the bandgap, only electrons at the edge of the valence band can interact with it, resulting in a low absorption. However, when photons have more energy, they can interact with more electrons, resulting in a larger absorption.

Absorption depth According to equation (2.1), one can see that the absorption coefficient depends on the wavelength. Consequently, light penetrates differently into a semiconductor depending on its wavelength, i.e. the light does not travel the same distance in the material before being absorbed. From the absorption coefficient, it is easy to determine the distance at which light intensity falls below 36% of its original one (or to be reduced by a factor 1/e, with e the elementary charge). This distance is the absorption depth and corresponds to α -1 . It reveals that short wavelength light (high energy light) is absorbed over a short distance because of its higher absorption coefficient.

Shockley-Queisser limit

The Shockley-Queisser limit sets an upper bound for single-pass solar cell efficiency. It was calculated by William Shockley and Hans-Joachim Queisser in 1961 [START_REF] Shockley | Detailed Balance Limit of Efficiency of pn Junction Solar Cells Additional information[END_REF]. This limit predicts the maximum photovoltaic conversion efficiency using a single bandgap semiconductor, which is around 33.7% for a single p-n junction cell under typical sunlight conditions. It is calculated by considering the amount of electrical energy extracted per incident photon. The general approach to compute the Shockley-Queisser limit is to calculate the absorption flux and the flux emitted from the solar cell. The difference between these two is the current from the solar cell.

Overview of solar cell technology

One way to maximize absorption of the whole solar spectrum is to stack layers of different materials, each absorbing its specific part of the spectrum. In this context, new types of solar cell designs have emerged.

Crystalline silicon (c -Si) is the most matured semiconductor technology used for the production of solar cells. It relies on a thick wafer of a few hundred micrometers in order to absorb sunlight. The state-of-art c -Si solar cell efficiency is 26.7% [START_REF] Yamamoto | High-efficiency heterojunction crystalline si solar cells[END_REF] Thin-film technologies only use few micrometers layers to absorb sunlight thanks to direct bandgap materials. III-V semiconductors are one of the best technologies to achieve high-efficiency thin-film solar cells. Thanks to the Shockley-Queisser limit, GaAs (Gallium arsenide) is the most efficient material for single-junction solar cells, with an efficiency of up to 28.8%. One way to maximize absorption of the whole solar spectrum is to stack different materials, each absorbing its part of the spectrum. The bottleneck for III-V cells is their high cost due to slow crystal growth and expensive III-V crystalline substrates. These materials are mainly used for spatial engineering.

Crystalline silicon modules have reduced the market price below 1 USD per Watt since 2013. This is due to improved efficiency and thinner c-Si wafers (nowadays about 160 ≈ 190 µm, in contrast to 300 µm in the early 2000s). In 2017, the worldwide PV module production was estimated at around 97.5 GWp. 95% of the total production is from wafer-based silicon technology (62% multi-Si and 33% mono-Si) and 5% of the total production is from thin-film technologies (CdTe, CIGS, and a-Si). III-V multijunction and concentrating photovoltaic (CPV) systems still account for very few percentage in the terrestrial photovoltaic application. Figure 2.6 displays the efficiency of research solar cells classified by materials and types from 1975 to nowadays.

The majority of solar cells suffer from losses, meaning that the absorption of sunlight is incomplete over the whole spectrum. This can be caused by the physical and electrical properties of the PV material. For example, Amorphous Silicon (a-Si) is a bad absorber in the red part of the spectrum (Λ ≈ 700 nm). This loss mainly comes from thermalization loss. Another common reason is the high reflection of the front surface of the solar cell.

Efficiency/cost trade-off

In order to increase the efficiency of PV cells, different strategies are possible. For example, one can increase the amount of absorber inside the cell by thickening the absorbing layers. This possibility is not privileged nowadays because it can considerably increase the cost of the cell. A more elegant way is to use nano-texturing and light trapping techniques, so that the light can be focused in the absorber layers as in [DTM + 17]. For such schemes, the standard theory developed by Eli Yablonovitch [START_REF] Yablonovitch | Statistical ray optics[END_REF], shows that the absorption enhancement factor has an upper limit of 4n 2 (knows as the Yablonovitch limit) where n is the refractive index of the absorber. For nanophotonic films with thickness comparable or even smaller than wavelength scale, ray optics picture, and some of the underlying assumptions in the conventional theory are no longer applicable. In that case, it can be shown that the absorption enhancement factor can go far beyond the Yablonovitch limit with proper design [MC12, BHS + 12]. In order to model the absorption of light in solar cells, one needs to solve the system of Maxwell's equations (section 2.2).

Maxwell's equations

The complete set of Maxwell's equations in their modern form, that describe the spatio-temporal evolution of electromagnetic waves on a domain Ω ⊂ R 3 , over a given time interval [0, T ] ,T > 0, are defined for all (x, t) ∈ Ω × [0, T ] by

∇ × E = - ∂B ∂t , (2.2)
∇ × H = ∂D ∂t + J, (2.3)
∇ • D = ρ, (2.4)
∇ • B = 0. (2.5)
Maxwell's equations with E and B the electric field and the magnetic induction, which are vectors in R 3 . Most of the physics textbooks considered E and B as "fundamental fields". We also introduce the electric displacement D, the magnetic field H, as well as the density of free electric charges ρ, and the free electric current density J. These quantities are all dependent on position x = t (x, y, z) and time t.

Constitutive relations

To close the system (2.2 -2.3), relations between (E, B) and (D, H) are required. In the general case, the constitutive relations are:

D = εE, B = μH, (2.6)
where ε and μ are tensors depending on x, t, E and B. One can make some assumptions on the considered material to temporarily simplify this system. If the material is linear, ε and μ are independent of E and B; if it is isotropic, ε ≡ εI 3 and μ ≡ µI 3 ; if it is homogeneous ε and µ are constant. Assuming all those assumptions, (2.6) becomes:

D = εE, B = µH.
It is standard to use relative permittivity ε r and relative permeability µ r , which are the normalization of ε and µ by ε 0 and µ 0 , respectively the vacuum permittivity and permeability. Of course, in vacuum we have ε r = 1 and µ r = 1. Hence, the constitutive relations are written as follows:

D = ε 0 ε r E, B = µ 0 µ r H. (2.7)
and Maxwell's equations for linear, homogeneous, isotropic, non-dispersive materials:

∇ × E = -µ 0 µ r ∂H ∂t , ∇ × H = ε 0 ε r ∂E ∂t + J.
(2.8)

Before being able to solve it, we need to add boundary and initial conditions to system (2.8). For numerical processing, we are adimensionning the system (2.8). The adimenssionated variables of (E, H, J) are noted ( E, H, J), and are computed as:

H = Z 0 H, E = E, t = c 0 t, J = Z 0 J,
with Z 0 = µ 0 ε 0 the vacuum impedance and c 0 = 1 √ ε 0 µ 0 the speed of light in vacuum. Hence, the normalized Maxwell's system is:

µ 0 c 0 Z 0 ∂ H ∂ t = - 1 µ r ∇ × E, ε 0 c 0 Z 0 ∂ E ∂ t = 1 ε r ∇ × H -J ,
Given the definitions of c 0 and Z 0 , one sees that µ 0 c 0 Z 0 = ε 0 c 0 Z 0 = 1. We can now drop the tilde notation to obtain the normalized Maxwell system, which will be used from now on.

∂H ∂t = - 1 µ r ∇ × E, (2.9)
∂E ∂t = 1 ε r (∇ × H -J) .
(2.10)

Maxwell normalized system

Dispersion models

Dispersion is a common phenomenon for all types of waves that travel through a medium. Wave's propagation is affected depending on the medium. For instance, for a polychromatic wave, all the wavelengths do not travel at the same speed through the medium. Among the numerous phenomena encountered in electromagnetics, many rely on the dispersive properties of materials. Dispersion is modeled by a frequency-dependent (wavelength-dependent) permittivity function ε(ω) (ε(λ)), often obtained from physical considerations. Regarding nanophotonics and photovoltaics applications, an accurate modeling of the ε for semiconductors in the visible spectrum is crucial. Several models aim to simulate this material specificity, as the Drude model [START_REF] Drude | Zur elektronentheorie der metalle[END_REF] and the Generalized dispersive model [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the fdtd method[END_REF][START_REF] Viquerat | Fitting experimental dispersion data with a simulated annealing method for nano-optics applications[END_REF].

Drude and Drude-Lorentz model

The Drude model is based on the kinetic theory of gases [START_REF] Drude | Zur elektronentheorie der metalle[END_REF]. This approximation consists of considering the metal as a static lattice of ions immersed in a free electrons gas. The interactions of these electrons with the ions are managed with collision frequency parameter γ d . With regard to electron-electron interactions, they are totally neglected. For the electron gas, this leads to the following classical equation of motion:

∂ 2 x ∂t 2 + γ d ∂x ∂t = - e m e E(t),
with m e the electron mass and e the electronic charge. Considering a harmonic time-dependence of the form e -iωt , one obtains:

x = e m e 1 ω 2 + iωγ d E.
By using the definition of the polarization P = -n e e x, with n e the electronic density, the latter equality can be rewritten as:

P = -ε 0 n e e 2 m e ε 0 1 ω 2 + iωγ d E,
and thus, the electric displacement becomes:

D = ε 0 1 - ω 2 d ω 2 + iωγ d E,
where

ω d = nee 2
meε 0 is called the plasma frequency of the electrons. It is common to include an additional parameter, χ b , describing the contribution of the bound electrons at infinite frequency [START_REF] Maier | Plasmonics: Fundamentals and Applications[END_REF]:

D = ε 0 1 + χ b - ω 2 d ω 2 + iωγ d E.
Then, the definition of the relative permittivity function is directly obtained by matching the previous expression with D = ε 0 ε r (ω) E:

ε r,d (ω) = ε ∞ - ω 2 d ω 2 + iωγ d , (2.11)
where ε ∞ = 1 + χ b is the permittivity at infinite frequency. Hence, in the same classical approach as before, a spring term is added to the equation of motion:

∂ 2 x ∂t 2 + γ l ∂x ∂t + ω 2 l x = - e m e E(t).
Following the same development as for the Drude model, one obtains the expression of a Lorentz pole:

ε r,l (ω) = - ∆εω 2 l ω 2 -ω 2 l + iωγ l .
The Drude-Lorentz model is then obtained by summing Drude and Lorentz terms:

ε r,dl (ω) = ε ∞ - ω 2 d ω 2 + iωγ d - ∆εω 2 l ω 2 -ω 2 l + iωγ l . (2.12)
with ∆ε a parameter representing the amplitude of the associated Lorentz pole.

Generalized dispersive model

Drude and Drude Lorentz models aims at modeling metals. For more specific materials, one usually cannot fit the dispersion profiles with the classical models. A more general approach is needed, allowing the fitting of any sort of optical measurements. A Padé style approximation is a simple analytical method for fitting an empirical set of points representing a material's permittivity function. The fundamental theorem of algebra allows to expand this approximation as a sum of a constant, one zero-order pole (ZOP), a set of first-order generalized poles (FOGP), and a set of second-order generalized poles (SOGP), as: Most of the standard dispersion models are naturally includes in this single formulation (Debye, Drude and Drude-Lorentz, Sellmeier's law). This model has been studied in [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF], and was also used by Vial et al in [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the fdtd method[END_REF] and by Han et al for the Complex-Conjugate Pole-Residue Pairs model (CCPRP) (see [START_REF] Han | Model dispersive media in fdtd method with complex-conjugate pole-residue pairs[END_REF]).

ε r,g (ω) = ε ∞ - σ iω - l∈L 1 a l iω -b l - l∈L 2 c l -iωd l ω 2 -e l + iωf l , ( 2
The strategy to fit the coefficients of (2.13) to experimental data is taken from [START_REF] Viquerat | Fitting experimental dispersion data with a simulated annealing method for nano-optics applications[END_REF]. This fitting is a not trivial optimization problem. Especially for an increasing number of poles (L1 and L2), one can be left with a significant maximization problem, which may have many local maxima. Simulated Annealing (SA) methods have proved to be particularly well suited and efficient for finding global maxima in these situations ( [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]). In practice, a set of experimental data (measured by physicists or taken from well-known database such as1 ) is provided to the SA algorithm, as well as a number of FOGP and SOGP. The spectrum of interest is a crucial factor in the quality of the fitting, and it should be chosen in accordance with the part of the light considered. Indeed, the number of poles needs to depend on the behavior of the experimental permittivity function in the selected frequency range. In this manuscript, the spectrum of interest is in general set to [300, 1000] nm, which constitutes a good representative range for PV applications. We show the fitting obtained for microcrystalline silicium (µc-Si) using 2 FOGP and 5 SOGP in figure 2.8. The optical measures are perfectly fitted over the whole range, with a relative error less than 1%. In chapter 4, a fit issue concerning bandgap modeling which is impacting the physical quantities specific to photovoltaics will be discussed.

Maxwell-GDM equations

In this section, the modified Maxwell's equations for dispersive media are derived by considering the previously introduced generalized dispersive model. One can derive the system of PDEs, following the steps described in [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF] for the Maxwell-Drude equations. The obtained system of PDEs, accounting for the generalized dispersive model (GDM) in time-domain is given by:

∂H ∂t

= -∇ × E, ∂E ∂t = 1 ε ∞   ∇ × H -J s -J 0 - l∈L 1 J l - l∈L 2 J l   , J 0 = (σ + l∈L 2 d l )E, J l = a l E -b l P l ∀l ∈ L 1 , ∂P l ∂t = J l ∀l ∈ L 1 , ∂J l ∂t = (c l -d l f l )E -f l J l -e l P l ∀l ∈ L 2 , ∂P l ∂t = d l E + J l ∀l ∈ L 2 .
(2.14)

Maxwell-Generalized dispersive model

We now switch to the numerical treatment of the previously stated systems of PDEs. The basic solver that we consider in this work combines a discontinuous Galerkin (DG) approach for the discretization in space (sections 2.4 and 2.3.3) with a fully explicit scheme for time integration scheme (section 2.5), yielding a so-called DGTD method.

Discontinous Galerkin formulation for Maxwell equations

In this section, we present the main principles and numerical ingredients of the DGTD method formulated on a conforming tetrahedral mesh and assuming an uniform interpolation order across all the elements of the mesh, which is essentially the outcome of the Ph.D. thesis of J. Viquerat [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF]. In the next chapters of this manuscript, we will extend this method and explain in more detail our methodological contributions motivated by the study of light interaction with nanostructured materials for applications relevant to photovoltaics. We consider the case of Maxwell equations in non-dispersive materials for the presentation of the different aspects of the DG formulation in space (sections 2.4.1 to 2.6, and then elaborate on the main modifications for dealing with dispersive materials.

Weak formulation

Let Ω ⊂ R 3 be a bounded convex domain, and n the unitary outward normal to its boundary ∂Ω.

Let Ω h be a discretization of Ω relying on a quasi-uniform triangulation T h (∃δ, ∀T i ∈ T h , ∀k ∈ V i , h k ≤ δh i , where h i is the size of the element T i ). This triangulation verifies T h = N i=1 T i , where N ∈ N * is the number of mesh elements and (T i ) i∈ 1,N is the set of elements of T h . We denote by a ik = T i T k the internal faces of between the adjacent cells T i and T k , and n ik is the unit normal vector to a ik , oriented from T i to T k . For each cell T i , we define as V i the set of indices k ∈ 1, N , such that T i and T k are sharing the face a ik .

From now on, we can write the weak formulation of problem (2.9) locally in the cell T i . Taking the L 2 scalar product of each term with a vector test function ψ, we obtain the following variational problem:

Find

(E, H) ∈ H 0 (curl, Ω h ) × H (curl, Ω h ) such that ∀ ψ ∈ H (curl, Ω h ) , ˆTi µ r ∂H ∂t • ψ + ˆTi ∇ × E • ψ = 0, ˆTi ε r ∂E ∂t • ψ - ˆTi ∇ × H • ψ = - ˆTi J • ψ.
Then, using classical vectorial calculus and Green formulas gives:

ˆTi µ r ∂H ∂t • ψ + ˆTi E • ∇ × ψ = ˆ∂T i (ψ × E) • n i , ˆTi ε r ∂E ∂t • ψ - ˆTi H • ∇ × ψ = - ˆTi J • ψ - ˆ∂T i (ψ × H) • n i .
One can notice that the previous equality only holds if the boundary terms ´∂T i (ψ × H) • n i exist. By using the properties of the mixed product, the previous system becomes:

(ψ × E) • n i = (E × n i ) • ψ.
This implies that taking E in H 0 (curl, Ω h ) requires the existence of the trace of ψ on ∂T i . Thus, we choose to take

ψ in H 1 (Ω h ). Therefore, ∀T i , ∀ψ ∈ H 1 (Ω h ), ˆTi µ r ∂H ∂t • ψ + ˆTi E • ∇ × ψ = ˆ∂T i (E × n i ) • ψ, ˆTi ε r ∂E ∂t • ψ - ˆTi H • ∇ × ψ = - ˆTi J • ψ - ˆ∂T i (H × n i ) • ψ.
(2.15)

Local weak formulation

Discretization in space

We start by defining the approximation space V h such that:

V h = v ∈ L 2 (Ω) 3 , v T i ∈ (P p (T i )) 3 ∀ T i ∈ T h ,
where P p (T i ) is the space of polynomial functions of maximum degree p on T i . We denote by (H h , E h , J h ) the semi-discrete fields seeked in space V h , and there restrictions on the cell T i as

(H i , E i , J i ) = H h T i , E h T i , J h T i .
We also define a set of scalar basis functions (φ ik ) 1 k d i for each T i , with d i the number of degrees of freedom (d.o.f.) per dimension. Additionally, we associate the three vectors φ v ik to each scalar basis function, such that:

φ 1 ik =   φ ik 0 0   , φ 2 ik =   0 φ ik 0   , φ 3 ik =   0 0 φ ik   .
For a 3D system, we notice that E i and H i are two vectors of three components:

E i =   E x i E y i E z i   , H i =   H x i H y i H z i   ,
which can both be locally expanded on the chosen set of basis functions:

E v i = d i j=1 E v ij φ ij , H v i = d i j=1 H v ij φ ij , v ∈ {x, y, z} . (2.16)
Therefore, for practical purpose, we define the three vectors of d i components:

E v i =    E v i1 . . . E v id i    , H v i =    H v i1 . . . H v id i    , v ∈ {x, y, z} ,
as well as the following 3d i components vector:

E i =       E x ij 1 j d i E y ij 1 j d i E z ij 1 j d i       , H i =       H x ij 1 j d i H y ij 1 j d i H z ij 1 j d i       .
Those vectors will be used to write the matrix-vector semi-discrete variational formulation of our system, analogously to (2.15). However, before that, the boundary terms in this formulation require some particular treatment before going further on in the spatial discretization process.

Numerical fluxes

Given that the test functions ψ can now be discontinuous at the interfaces a il (between cells T i and T l ), it is important to notice that the surface integrals, define as:

ˆail (E h × n il ) • ψ, and
ˆail (H h × n il ) • ψ, (2.17)
are not clearly define. Indeed, electromagnetic fields E h and H h can relate to either the field value in T i or in T l . We need to define the numerical fluxes (or traces), which will allow to recover a proper definition of the surface integrals. This will permit to reconnect the field values between neighbouring cells.

Remark : We specify that the numerical fluxes do not have a unique valid choice. In fact, different choices can lead to stable and convergent discrete schemes, with different numerical properties.

We replace in equations (2.17) E h and H h by E * and H * such that:

ˆail (E * × n) • ψ, and
ˆail (H * × n) • ψ.
It is well known that the flux can be seen as the solution of a Riemann problem at cell interfaces.

In this thesis, we will exploit two very common flux choices. The first one is the centered flux, which is defined as:

E * = E i + E l 2
, and

H * = H i + H l 2 . (2.18)
This flux is non-dissipative and leads to an L 2 spatial convergence in h p . When it is coupled to a non-dissipative time-integration scheme (for example Leap-Frog), this yields to a non-dissipative DGTD scheme [START_REF] Fezoui | Convergence and stability of a discontinuous galerkin time-domain method for the 3d heterogeneous maxwell equation on unstructured meshes[END_REF]. The second numerical flux considered is the upwind flux, which is defined as:

E * = 1 Y i + Y l ({Y E} il + αn × H il ) , and
H * = 1 Z i + Z l ({ZH} il -αn × E il ) , (2.19)
where {A} il = A i + A l is twice the mean value of A at the interface, A il = A l -A i is the jump of A at the interface, and α ∈ [0, 1] is a tunable parameter that allows to vary between the centered flux (2.18) taking α = 0, to a fully upwind flux for α = 1. The jump term of the upwind flux introduces dissipation in the DG scheme. Those dissipation can help damping nonphysical modes when instabilities occur [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF]. It also leads to a better L 2 spatial convergence (as h p+1 against h p for the centered flux).

DG matrices

In this part, we present the FE matrices which will allow us to write the matrix-vector form of system 2.14. We choose the test functions ψ such that they are the 3 d i vectors φ v ik , which constitutes the Galerkin choice:

ˆTi µ r ∂H i ∂t • φ v ik + ˆTi E i • ∇ × φ v ik = l∈V i ˆail (E * × n il ) • φ v ik , ˆTi ε r ∂E i ∂t • φ v ik - ˆTi H i • ∇ × φ v ik = - l∈V i ˆail (H * × n il ) • φ v ik - ˆTi J i • φ v ik .
(2.20)

Mass matrix We start by considering the time-derivative term of the E evolutionary equation of system (2.20). The x component is:

ˆTi ε r ∂E i ∂t • φ x ik = ˆTi ε r ∂E x i ∂t φ ik = d i j=1 ∂ ∂t E x ij ˆTi ε r φ ij φ ik = M εr i ∂ E x i ∂t k
where M εr i is the mass matrix, of dimension d i × d i :

(M εr i ) jk = ˆTi ε r φ ij φ ik , with (j, k) ∈ 1, d i 2 .
Stiffness matrices We now focus on the curl integral of the E evolutionary equation of system 2.20. The x component is:

ˆTi H i • ∇ × φ x ik = ˆTi H y i ∂φ ik ∂z -H z i ∂φ ik ∂y = ˆTi d i j=1 H y ij φ ij ∂φ ik ∂z -H z ij φ ij ∂φ ik ∂y = d i j=1 H y ij ˆTi φ ij ∂φ ik ∂z - d i j=1 H z ij ˆTi φ ij ∂φ ik ∂y = K z i H y i -K y i H z i k = -K i × H i x k .
with the three stiffness matrices defined as:

(K v i ) jk = ˆTi φ ij ∂φ ik ∂v for v ∈ {x, y, z} , with (j, k) ∈ 1, d i 2 .
We can also define the global 3d i × d i stiffness matrix that will be used in the matrix-vector form of the final system:

K i =   K x i K y i K z i   ,
Flux matrix Finally, we consider the right member of system 2.20 that contains the flux contribution. Here we are using the centered flux (2.18), but the upwind case, as well as the generalization to other fluxes is straightforward. We also note that in the conforming case, the field expansion defined on a il over the basis functions of cells T i or T l is equivalent. We proceed as previously, focusing on the x component of the flux of the E evolutionary equation of system 2.20:

ˆail (H * × n il ) • φ x ik = ˆail H y * n z il -H z * n y il φ ik = ˆail H y i + H y l 2 n z il - H z i + H z l 2 n y il φ ik = 1 2 d i j {H y } il n z il -{H z } il n y il ˆail φ ij φ ik = S il H * × n il x k
where the flux matrices are:

(S il ) jk = ˆail φ ij φ ik , with (j, k) ∈ 1, d i 2 .
Remark : In non-conforming cases (which includes hybrid meshes, non-conforming meshes and local approximation order), one should consider the proper expansions of the field over both cells respective bases. This is presented in details in the chapter 3.

Semi-discrete formulation Thanks to the definition of the FE elementary matrices, we can define both global mass and flux matrices, which will allow to write the semi-discrete formulation in a compact form:

M u i =   M u i 0 d i ×d i 0 d i ×d i 0 d i ×d i M u i 0 d i ×d i 0 d i ×d i 0 d i ×d i M u i   , S il =   S il 0 d i ×d i 0 d i ×d i 0 d i ×d i S il 0 d i ×d i 0 d i ×d i 0 d i ×d i S il ,   .
which leads to the following matrix vector expression of the semi-discrete DG scheme for Maxwell's equations:

M µr i ∂ H i ∂t = -K i × E i + l∈V i S il E * × n il , M εr i ∂ E i ∂t = K i × H i - l∈V i S il H * × n il -M i J i .
(2.21)

Semi-discrete scheme

Elements mapping

One of the strength of the DG method is that the FE matrices described in 2.4.4 are not stored for each element of T h , but are calculated only once for all on a reference element noted T and then mapped on the considered physical tetrahedron T i . Let T be defined as follows in the ξ = (ξ, η, ζ) coordinates system:

T = (ξ, η, ζ) ∈ R 3 + , ξ + η + ζ ≤ 1 .
Then we define the physical tetrahedron in the x = (x, y, z) coordinates system as the image of T by a mapping ψ T i (see figure 2.9) defined as:

ψ T i : T → T i , such that, ∀ξ ∈ T , x = ψ T i (ξ) .
The vertices of T are noted (A 1 , A 2 , A 3 , A 4 ), and the vertices of T i are noted (v 1 , v 2 , v 3 , v 4 ). In this case, ψ is a linear combination of ξ, η and ζ defined as:

ψ T i (ξ) = v 1 + (v 2 -v 1 )ξ + (v 3 -v 1 )η + (v 4 -v 1 )ζ.
Let (φ ij ) j=1..d i be the basis functions on T i , and φ j j=1..d i defined by we can write the mass matrix on the element T i as:

φ j = φ ij • ψ T i on T . Then A 1 A 2 A 3 A 4 ξ ζ η v 1 v 2 v 3 v 4 x z y x = ψ T i (ξ)
(M i ) jk = ˆTi φ ij (x)φ ik (x)dx = ˆ T φ j (ξ) φ k (ξ) J ψ T i dξ,
where J ψ T i (ξ) is the jacobian matrix of the mapping ψ, defined as:

J ψ T i jl = ∂x j ∂ξ l jl =   (v 2 -v 1 ) x (v 3 -v 1 ) x (v 4 -v 1 ) x (v 2 -v 1 ) y (v 3 -v 1 ) y (v 4 -v 1 ) y (v 2 -v 1 ) z (v 3 -v 1 ) z (v 4 -v 1 ) z   .
In this case, the determinant J ψ T i is constant and only depend on (v 1 , v 2 , v 3 , v 4 ). Hence, the mass matrix for each physical tetrahedron of T h is a simple multiplication of the mass matrix calculated on the reference tetrahedron:

(M i ) jk = J ψ T i M jk .
Similar situation occurs for stiffness and flux matrices by using the following change of variables (see [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF] for additional details):

(K v i ) jk = ˆTi φ ij (x) ∂φ ik ∂v (x) dx = ˆ T φ j (ξ) J ψ T i J -1 ψ T i ∇ ξ φ k (ξ) v dξ = J ψ T i J -1 ψ T i ˆ T φ j (ξ)∇ ξ φ k (ξ)dξ v = 3 m=1 J ψ T i J -1 ψ T i vm ˆ T φ j ∂ φ k ∂ξ m = 3 m=1 J ψ T i J -1 ψ T i vm K m jk ,
where ∇ ξ is the gradient operator in the T basis. Therefore, the stiffness matrix on T i can be built from the stiffness matrices K m jk calculated on the reference element. Finally, the flux matrix can be rewritten as:

(S il ) jk = ˆail φ ij φ lk ds = ˆ t φ j φ k J ψ T i J -1 ψ T i n d s = J ψ T i J -1 ψ T i n S jk ,
with t the reference triangle. As for the other matrices, these matrices can be deduced from the S jk . The unitary normals also transform as:

n il = J -1 ψ T i n J -1 ψ T i n .
Therefore, even for the flux matrix, we only need to store a single (3 × 3) matrix instead of once per cell.

Remark : If one is using local polynomial orders, flux matrices must be computed for each possible order's jump. If the mesh is non-conform or hybrid, flux matrices must be computed for all the non-conformal interfaces.

Polynomial expansion basis

In this part we discuss the choice of the basis function φ j j=1..d i . Over all the possible polynomial bases available, Lagrange polynomials are quite a common choice. They can be defined by a set of interpolation nodes distributed across the cell. The Lagrange interpolants L i are defined by the following property:

L i (x j ) = δ ij , ∀(i, j) ∈ 1, d i 2 .
There must be the same number of nodes and polynomials in order to define a complete basis. In a tetrahedron the number of Lagrande nodes needed to have a polynomial order p is equal to:

n(p) = (p + 1)(p + 2)(p + 3) 6 ,
and for each of its faces, it is:

s(p) = (p + 1)(p + 2) 2 .
We choose to use Lagrange polynomials with equispaced nodes as it allows a simple integration of the FE matrices on the reference element since the nodes positions are known exactly.

Remark : In the case of an hexahedral cell, the number of Lagrange nodes requires for a polynomial order q is: n(q) = (q + 1) 3 and s(q) = (q + 1) 2 .

DG operator spectrum

With DG, the global matrix of the spatial operator can be assemble such as FEM. Compared to the FEM global matrix, the DG one is of much larger dimension, even for the same number of cells, and the same interpolation order. This is due to the degree of freedom (d.o.f.) which are duplicated on the cell's interfaces. However, unlike the FEM global matrix, local matrices do not overlap. Indeed, the degrees of freedom exist independently at each cell interfaces on both sides. Eventually, we can rewrite the semi-discrete DG problem as the following generic form:

∂W h ∂t (t) = A h W h (t) + B h (t), (2.22)
where B h (t) represents the exterior sources, and the matrix A h holds all the information relative to the spatial discretization.

Time discretization

In section 2.4, the semi-discrete scheme was obtained (2.21) by discretizing the spatial derivatives of Maxwell's equations. Similarly, the time derivatives need to be considered for the time discretization. Then, we consider the reduced problem of the 1D Maxwell's equations:

µ r ∂H y ∂t = ∂E z ∂x , ε r ∂E z ∂t = ∂H y ∂x + j(t).
By dropping the spatial subscripts of the unknowns, the semi-discrete formulation associated to this system becomes:

M µr i ∂H i ∂t = K i E i + [E * ] x i x i-1 , M εr i ∂E i ∂t = K i H i + [H * ] x i x i-1 -M i J i , (2.23)
with:

E * = 1 Y i + Y l ({Y E} il + α H il ) , H * = 1 Z i + Z l ({ZH} il + α E il ) .
In the previous equalities, subscript l designates i -1 or i + 1. When trying to go from system (2.23) to system (2.22), a non-symmetric A matrix is obtained in the general case:

A = A α,H A H A E A α,E .
The off-diagonal blocks represent both the centered part of the flow and the stiffness part, while the diagonal blocks represent the upwind contribution. Therefore, A is strictly anti-diagonal in the case of centered flows. Furthermore, the following properties are obtained in the case of a homogeneous medium To begin with, we consider the simple case of vacuum (ε r = µ r = 1). In this case, one has a symmetric A matrix and thus a diagonalizable system. Thus, one can just retain the corresponding formulation in the diagonalized basis, for the study of the time-stepping schemes. This has the effect of reducing to a system of ODEs of the form:

(Y i = Y l = Y ): 1. A α,
∂φ ∂t (t) = λφ(t) + b(t) ≡ f (t, φ(t)), (2.24)
for each λ eigenvalue of A. A number of time-integration techniques from the ODE community are appropriate to solve (2.24). We can classified time-stepping methods into two main categories: the explicit time integration techniques, and the implicit methods. For the first category, the time state φ(t + ∆t) is computed explicitly from φ(t). For the second category, the time-updated solution is obtained by solving an implicit expression of the form g(φ(t), φ(t + ∆t)) = 0. It leads to the resolution of a linear system of equations at each time-step. As expected, the implicit method is more expensive than explicit technique. But, implicit methods are generally unconditionally stable. This means that any choice of ∆t leads to a stable algorithm. For explicit methods, a numerical criterion on the time-step called the Courant-Friedrichs-Lewy (CFL) condition, must be respected. Otherwise, the resulting algorithm will be unstable and blow up. In conclusion, explicit time-stepping generaly requires much more time-steps but each time-step is significantly less numerically costly. Hereafter, solutions are investigated on intervals of the form [0, T ] with T > 0, discretized in time-steps of length ∆t. In order to refer the discrete approximation of φ(t n ), with t n = n∆t, we used the notation φ n . A family of DG methods, called space-time DG methods [PFT00] [START_REF] Angulo | An explicit nodal spacetime discontinuous galerkin method for maxwell's equations[END_REF], is designed to deal with time derivative similarly to space. The main disadvantage of these methods is that they usually lead to an implicit scheme. In the next part, we present the Runge-Kutta (RK) and Low-Storage Runge-Kutta (LSRK) schemes.

Classical RK schemes

Runge-Kutta time schemes are a class of multi-stage algorithms based on the multiple evaluations of the right side of (2.24) to evolve the system in time. Suppose that problem (2.24) is integrated between t and t + ∆t:

φ(t + ∆t) = φ(t) + ˆt+∆t t f (u, φ(u))du.
(2.25)

We can approximate (2.25) using a quadrature formula with s terms:

φ(t + ∆t) φ(t) + ∆t s j=1 β j f (t + δ j ∆t, φ(t + δ j ∆t))du, (2.26)
where (β j ) j∈ 1,s and (δ j ) j∈ 1,s are constants depending of the choice of the quadrature formula.

In order to evaluate the φ(t + δ j ∆t) values, RK methods use a so-called prediction/correction technique. It builds on the previous guesses to calculate the next one. The n th time-step with an s-stage RK algorithm is can be written as follow:

φ 1 = f (t n , φ n ) φ k = f   t n + δ k ∆t, φ n + ∆t s j=1 α j,k φ j   for k = 2, ..., s, φ n+1 = φ n + ∆t s j=1 β j φ j , (2.27)
where we suppose that φ 0 = φ n . The system (2.27) is implicit in the general case. The summation in each intermediate stage extends to the maximum number of stages. For practical and efficacity reason, we choose to work with explicit RK schemes (i.e. α j,k = 0 ∀k ≥ j). They can be written as:

φ 1 = f (t n , φ n ) φ k = f   t n + δ k ∆t, φ n + ∆t k-1 j=1 α j,k φ j   for k = 2, ..., s, φ n+1 = φ n + ∆t s j=1 β j φ j .
(2.28)

Explicit RK schemes

For larger number of stages, the principle of a RK algorithm remains similar to the two stages algorithm. At time-step n, the first stage of the method is:

φ 1 = f (t n , φ n ) ≡ f tn , (2.29)
and the second stage is:

φ 2 = f (t n + δ 2 ∆t, φ n + ∆tα 1,2 φ 1 ) . (2.30)
In (2.30), one can recognize φ n + ∆tα 1,2 φ 1 as a forward Euler method, that provides a firstorder approximation of φ(t n + α 1,2 ∆t) by exploiting the approximate value of the slope at t = t n calculated during the stage (2.29). One can obtained φ 2 by setting δ 2 = α 1,2 . The method is done in two stages then, the solution at the next time-step is obtained as:

φ n+1 = φ n + ∆t (β 1 φ 1 + β 2 φ 2 ) = φ n + ∆t (β 1 f tn + β 2 f (t n + δ 2 ∆t, φ n + ∆tα 1,2 φ 1 )) = φ n + ∆t (β 1 f tn + β 2 f (t n + α 1,2 ∆t, φ n + ∆tα 1,2 φ 1 )
) . Using Taylor expansion to the first order in ∆t on the last term gives:

φ n+1 = φ n + ∆t (β 1 + β 2 ) f tn + β 2 α 1,2 ∆t 2 ∂f ∂t (t n , φ n ) + f tn ∂f ∂φ (t n , φ n ) + O ∆t 3 ,
that can be matched with the Taylor expansion of φ, which gives:

φ(t + ∆t) = φ(t) + ∆tf t + ∆t 2 2 ∂f ∂t (t, φ) + f t ∂f ∂φ (t, φ) + O ∆t 3 ,
where f t = f (t, φ(t)). In practice we drop the residual O ∆t 3 , that leads to the following system:

β 1 + β 2 = 1, β 2 α 1,2 = 1 2 .
(2.31)

One can see RK schemes as a quadrature rule for which the sum of the weights is one, as recall in the first equation of (2.31). One also remark that (2.31) is underdetermined, which means that there is not a unique set of coefficients verifying (2.31). A common choice to satisfy the latter system is:

α 1,2 = 1,
β 1 = β 2 = 1 2 . (2.32)
This choice leads to a second-order RK method:

φ 1 = f (t n , φ n ) φ 2 = f (t n + ∆t, φ n + ∆tφ 1 ) φ n+1 = φ n + ∆t 2 (φ 1 + φ 2 ) .
(2.33)

Here we summarize the algorithm which is illustrated in figure 2.10. An estimate of ∂φ ∂t (t n ) is given by f tn during the first step. Then, via a forward Euler step, we obtain a first-order approximation φ n+1 * = φ n + ∆tφ 1 of φ(t + ∆t). In the second step, φ n+1 * is used to estimate ∂φ ∂t (t n+1 ). Finally, a second forward Euler step is used to obtain φ n+1 .

Similarly, one can obtained the four-step fourth-order RK algorithm by exploiting the average of four different evaluation of the slope (figure 2.11). The most classical version of this scheme is:

φ 1 = f (t n , φ n ) φ 2 = f t n + ∆t 2 , φ n + ∆t 2 φ 1 φ 3 = f t n + ∆t 2 , φ n + ∆t 2 φ 2 φ 4 = f (t n + ∆t, φ n + ∆tφ 3 ) φ n+1 = φ n + ∆t 6 (φ 1 + 2φ 2 + 2φ 3 + φ 4 ) .
(2.34)

Remark : Obtaining an RK scheme of order p when p > 4 requires more than p stages ( [START_REF] Butcher | The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods[END_REF][START_REF] Lambert | Numerical Methods for Ordinary Differential Systems[END_REF]).

Low-storage RK schemes

Standard RK schemes can achieve high order accuracy. However, these algorithms needs to store all the s slope estimations in addition to the usual solutions. If this is not a major constraint for small problems and/or for a small number of substeps, this storage become a huge limitation for larger problems. Thus, a specific class of RK algorithms have been developed to overcome this issue. These algorithms require only two memory registers at every moment. They are called low-storage RK schemes, and they are usually presented under the Williamson formulation [START_REF] Williamson | Low-storage runge-kutta scheme[END_REF]:

φ 1 = φ n φ 2 = a k φ 2 + ∆tf (t n + d k ∆t, φ 1 ) , for k = 1, ..., s φ 1 = φ 1 + b k φ 2 , φ n+1 = φ 1 .
(2.35)

Low-storage RK schemes

Remark : The convergence to order p is not anymore ensured in s substeps. At least 5 stages are necessary to achieve a fourth-order convergence as see in [START_REF] Carpenter | Fourth-order 2n-storage runge-kutta schemes[END_REF].

The principle of LSRK schemes is slightly different from RK schemes. The final step is mixed within the slope estimation steps: starting from the initial point, the slope is approximated (with φ 2 ), and an Euler step is executed, giving an intermediate value of the solution (φ 1 ) at time station t n +b 1 . Then during the second stage, the slope stored in φ 2 receive a weighted average between the previous estimate and the new one (calculated from φ 1 at t n + b 1 ). This is illustrated in figure 2.12.

Time-step and CFL condition

The main drawback of explicit methods is the CFL condition. This time-step restriction is computed from an energy-based stability study. In this thesis, theoretical results from [START_REF] Fezoui | Convergence and stability of a discontinuous galerkin time-domain method for the 3d heterogeneous maxwell equation on unstructured meshes[END_REF] are used. Consequently, for a space discretization with polynomial order p, the time-step is chosen as follows:

∆t p = c p min T i ∈T h V T i A T i , (2.36)
where V T i is the volume, A T i is the area of cell T i , and c p is an order-dependent constant. The maximal acceptable value for c p , can be determined on a basic test case.

Domain truncation

Boundary conditions

Each computational setup needs to be closed by boundary conditions, which depends on the physics of the problem. Most of the time, photovoltaic problems are open. The considered physical domain must be artificially truncate due to the fact that computational domain cannot be infinite. In PV applications, some problems are peridiodic. This means they can be truncated using the so-called periodic boundary conditions (PBC). PBCs describe an infinite repetition of the same pattern. Periodicity can also be represented through symmetric domain, using perfect electric conductor condition (PEC) and perfect magnetic conductor condition (PMC).

In DGTD methods, the boundary conditions are imposed by adding a layer of cells called 'ghost cells'. This layer is outside the computational domain (see figure 2.13). Then, specific values of the fields are set in these ghost cells. Thus, the behavior of the solution on the boundary is controlled, without any special treatment, via numerical fluxes. Hereafter, the field inside the ghost cells are noted E gc and H gc , while the fields in the boundary cells are denoted E bc and H bc . Perfect electric conductor condition Perfect electric conductor condition (PEC) is an idealized material which exhibits infinite electrical conductivity. A zero tangential electric field is enforced by setting the field values in the ghost cells on a PEC boundary as:

(E gc , H gc ) Ω (E bc , H bc )
E gc = -E bc and H gc = H bc .
Perfect magnetic conductor condition Perfect magnetic conductor condition (PMC) reprensents the reciprocal of the PEC condition. Most of the time, it is used to impose symmetry planes.

A zero tangential magnetic field is enforced by setting the field values in the ghost cells on a PMC boundary as:

E gc = E bc and H gc = -H bc .
Absorbing boundary condition Absorbing boundary conditions (ABC) allow to partially absorb fields radiating out from the physical domain. It exists several forms of these conditions. The first-order Silver-Müller boundary condition is often used [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF]:

n × (E + Z (n × H)) = 0 n × (H -Y (n × E)) = 0 (2.37)

Silver-Müller boundary condition

Imposing the Silver-Müller boundary condition is similar to set the incoming flux to zero on the boundary. For this reason, its expression depends on the upwinding factor α:

Y gc = Y bc , Z gc = Z bc , E gc = 0 + 1 -α Y bc n × H bc and H gc = 0 - 1 -α Z bc n × E bc .
Remark : The Silver-Müller condition perfectly absorbs normally-incident plane waves, when imposed on the boundary. However, when waves are incident at increasing angles, its performance decreases.

Periodic boundary condition Periodic boundary conditions (PBC) can be used in order to simulate infinite mono-or bi-directional arrays, considering one elementary pattern. Cells from a periodic boundary face are matched with cells on the opposite boundary face of the domain. This way, every cell has a neighbor which is well-defined, and standard fluxes can be applied.

Remark : A periodic mesh is a prerequisite to use PBC. In other words, opposite faces in the periodic direction must match.

Perfectly matched layers

A novel numerical concept has been developed in 1994 by Bérenger in order to overcome the limitations of ABC. The objectives is to absorb the waves radiated from a system. To do so, Bérenger defined an artificial volume surrounding the physical domain in which a damping should occur progressively. Those artificial volumes are known as perfectly matched layers (PML). Figure 2.14 illustrates the functioning of PMLs. Outgoing wave propagates in the physical domain toward the PML, and crosses the interface. No reflection occurs, and the wave is progressively damped by the artificial medium while it continues to propagate in the PML. At one point, the wave will encounters the boundary of the computational domain, which is usually PEC or ABC. In both cases, the remaining of the wave will be totally (PEC) or partially (ABC) reflected toward the domain. In practice, the wave will travel a second time over the PML length, causing more damping. In general, when it enters again in the physical domain, the amplitude of the wave is attenuated by several orders of magnitude. For this reason, the error induced by PMLs is supposed PML Physical domain PEC No reflection occurs, and the wave is progressively damped by the artificial medium and continues to propagate in the PML. At one point, the wave might encounters the edge of the computational domain, which is usually PEC or ABC. In both cases, the remaining of the wave will be totally (PEC) or partially (ABC) reflected toward the domain. In practice, the wave will travel a second time over the PML length, causing more damping. In general, when it enters again in the physical domain, the amplitude of the wave is attenuated by several orders of magnitude.

to be small enough not to lose the benefits of high-order methods. ABCs can be seen as "geometric" condition (i.e. it only becomes more efficient with a larger distance from the source), while PMLs take advantage of the high-order spatial discretization to allow higher levels of damping.

PMLs have evolved since Bérenger's implementation, and several varieties are now available [START_REF] Bérenger | Perfectly matched layer for computational electromagnetics[END_REF]. One can cite the uniaxial PML (UPML) [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF], and the complex frequency-shifted PML (CFS-PML), which are in use for Maxwell's equations. The last one is used in this work.

CFS-PML

The complex frequency-shifted PML (CFS-PML) are based on a complex stretch of coordinates in the spatial operator. This approach was introduced by Kuzuoglu and Mittra [START_REF] Kuzuoglu | Frequency dependence of the constitutive parameters of causal perfectly matched layers[END_REF]. However further developments showed that this method was found successful for wrong reasons [START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF]. That being said, since [START_REF] Roden | Convolutional pml (cpml): An efficient fdtd implementation of the cfs-pml for arbitrary media[END_REF] the method has been broadly used in the FDTD community. In 2011, a DGTD implementation of this PML has been proposed by König [Kön11]. The next paragraph is based on this implementation. The following change of variables imposes the complex stretching:

∂ ∂x → 1 s x (ω) ∂ ∂x , ∂ ∂y → 1 s y (ω) ∂ ∂y , ∂ ∂z → 1 s z (ω) ∂ ∂z ,
with:

s k (ω) = κ k - σ k iω -α k , with k ∈ {x, y, z} ,
where σ k is the loss rate of the PML. Two new parameters are introduced: κ k that represents a real stretching factor, which effect is only to artificially lengthen the PML; α k which is the actual frequency shift, since it moves the pole of s k (ω) from ω p = 0 to ω p = iω. Linearly-growing instabilities in long-time computations can arise, when

α k = 0 [BPG02].

Sources

One of the major concern to simulate realistic configuration is a good control of the properties of incident fields. Indeed, the physical response of a nanophotonic system depends mainly on how it is excited. In this work, only plane waves are considered since they can model the solar light, even if several sources are available from realistic physical applications as laser beams [START_REF] Novotny | Principles of nano-optics[END_REF], or dipoles.

Plane waves

The most simple kind of source are plane waves. Indeed, they are commonly used in numerical electromagnetics to determine the fundamental properties of a physical system even if these waves correspond to an asymptotic physical configuration (i.e. any radiating source propagating on a sufficiently large distance should look like a plane wave). Resonances and modes of a physical system can be excited depending on the spectral profile of the source. Monochromatic plane wave is the most basic kind of time dependence. In the temporal domain, its expression is:

E(t) = E 0 sin(ω 0 t). (2.38)
One can notice that the Fourier transform of (2.38) is proportional to the Dirac function δ ω 0 . This type of source can be useful in several cases. However, running a whole time-domain simulation only to obtain the response of the system at one frequency may not be worth it, and a frequencydomain methods may be more suited for this type of source. In order to obtain a broader frequency spectrum, one can use:

E(t) = E 0 sin(ω 0 (t -t 0 )) e -(t-t 0 ) 2 2σ 2 , (2.39)
which is a Gaussian function of width σ centered around t 0 , and modulated by a sine function. Its Fourier transform can be computed as follow:

E(ω) = E 0 iσ π 2 e iωt 0 e -σ 2 (ω-ω 0 ) 2 2 -e -σ 2 (ω+ω 0 ) 2 2 , (2.40)
This means that such a pulse traveling through a structure will excite it on a wideband of wavelength, and not just on a single one. In our case, the idea is to use a pulse which contains the whole visible spectrum.

TF/SF formulation

Different possibilities are available in order to impose the plane waves inside the physical domain.

A first simple one is to use the ghost cells on the ABCs. This solution is no longer viable when PMLs are used (imposing fields in the PMLs will directly damp them). A possible alternative is to define an additional artificial contour inside the physical domain and between PMLs and the scatterer, on which the field could be imposed directly. For periodic domains, the artificial contour can be replaced with infinite artificial surfaces. First we consider the splitting of the electric field in two parts: where E tot is the total field, E inc is the incident field, and E sca is the scattered field. In our case, the incident field is a plane wave and it is known since it is imposed analytically. We consider now a splitting of the computational domain in two regions: the first one which in which the total field is computed, while in the second region the scattered field is computed. We call the interface between these two regions the total field/scattered field (TF/SF) interface. The DGTD formulation 2.21 is still valid in both region, and no modification is required, except for the computation of the fluxes between both regions. Consider a TF/SF interface between two cells, such that the local cell T i is in the total field region and its neighbor cell T l is the scattered field region. The upwind flux calculated for cell T i is:

E tot (x, t) = E inc (x, t) + E sca (x, t), (2.41) 0 2 4 6 8 10 -1 -0.5 0 0.5 1 t (a.u.) E(t) σ = 1 σ = 1.5 (a) Time-domain 0 1 2 3 4 5 6 0 10 20 30 ω E (ω) σ = 1 σ = 1.5 (b) Frequency-domain
E * ,tot = 1 Y i + Y l ({Y E tot } il + αn × H tot il ) , with {Y E tot } il = Y i E i,tot + Y l E l,tot and H tot il = H l,tot -H i,tot
. However, the field values corresponding to cell T l are not E l,tot and H l,tot , but E l,sca and H l,sca . Hence, the flux formulation must be modified to account for this difference. By considering (2.41), one easily shows that the right flux can be calculated as follows:

E * ,tot = 1 Y i + Y l ({Y E} il + αn × H il ) E * + 1 Y i + Y l (Y l E inc + αn × H inc ) E * ,inc , (2.42)
Symmetrically, the upwind flux for cell T l in the scattered field region is:

E * ,sca = E * -E * ,inc . (2.

43)

Remark : The fluxes computation at the TF/SF interface can be handled by a separate additional loop on the TF/SF faces during the fields update. causing a minimal overhead.

The TF/SF decomposition are interesting for several reasons. In the presence of PMLs, they can be used as an "interior boundary" to impose incident fields. It is possible to observe both the total field in near of the scatterer and the scattered field far from it in a single computation, and it helps to compute relevant scattering quantities, such as reflection, transmission and absorption (see chapter 4). In the previous chapter, we have presented a DG formulation using a uniform polynomial order across the whole computational domain, coupled with an explicit time-stepping scheme. However, when dealing with realistic and complex geometries, the resulting mesh can show significant variations in cell size. In such cases, the time-step is imposed by the smallest cells and depends on the polynomial order chosen, which can lead to prohibitive computational costs for high polynomial orders. Indeed, in such a scenario, high-order interpolation is only needed in the coarser cells, while both the computational overhead and the smaller time-steps in the finer cells have negative impacts.

In order to overcome this limitation, several possibilities can be considered. For example, one can replace the explicit time scheme by an implicit time-stepping algorithm. This results in the resolution of a large linear system at each (larger) time-step [START_REF] Verwer | Unconditionaly stable integration of maxwell's equations[END_REF]. Some hybrid implementations, called implicit-explicit (IMEX), were also proposed. For instance, in [START_REF] Verwer | Component splitting for semi-discrete maxwell's equations[END_REF] and [START_REF] Moya | Locally implicit discontinuous Galerkin time-domain methods for electromagnetic wave propagation in biological tissues[END_REF], the authors consider a hybrid formulation where the smallest cells are treated with an implicit scheme while using an explicit time integration for the bigger cells, thus limiting the time-step constraint. However, maintaining high-order time integration has proven to be difficult [START_REF] Moya | Locally implicit discontinuous Galerkin time-domain methods for electromagnetic wave propagation in biological tissues[END_REF]. Another possibility is to exploit local time-stepping (LTS). One of the main ideas is to divide the mesh into regions, which share the same time-step (see [START_REF] Piperno | Symplectic local time-stepping in non-dissipative dgtd methods applied to wave propagation problems[END_REF] and [START_REF] Diaz | Energy conserving explicit local time stepping for second-order wave equations[END_REF]). Moreover, high-order convergence is preserved by this method as shown in [START_REF] Diaz | Energy conserving explicit local time stepping for second-order wave equations[END_REF].

A complementary strategy, which is naturally suitable to DG methods, relies on the use of local polynomial orders. The above-mentioned drawbacks can be managed by imposing low orders in the smaller cells, and high-orders in the coarser cells. Doing so significantly reduces both the global number of degrees of freedom and the time-step restriction with a minimal impact on the accuracy of the method. Strategies exploiting locally adaptive (LA) formulations usually combine both h-and p-adaptivity in order to concentrate the computational effort in the areas of high field variations. In this chapter, we want to show that, starting from a given mesh with a homogeneous P k polynomial order, it is possible to exploit the LA strategy with all the polynomial orders P p (p ≤ k) and obtain a solution of similar accuracy for a reduced computational cost.

PV cells are made of an assembly of layers of materials. The difficulty in discretizing this type of domain lies in the interface between these layers, where highly detailed nanotextures are present. This area, which represents around 20 to 30% of the computational domain, is the main argument in favor of the use of an unstructured mesh. A strategy has already been proposed to drastically reduce the number of degrees of freedom as well as the computing time, notably for electromagnetism in [START_REF] Durochat | Méthode de type galerkin discontinu en maillages multi-éléments pour la résolution numérique des équations de maxwell instationnaires[END_REF], and more precisely for nanophotonics in [LVD + 14]. This strategy consists of using a structured hexahedral mesh in the homogeneous parts of the domain while using an unstructured mesh for the proper representation of geometry in the areas of importance.

In this chapter, we present two numerical strategies to deal with multiscale PV structures. First, a p-adaptive version of the DGTD method is studied in section 3.1. This is an extension of the one presented in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF]. Several strategies for the distribution of the interpolation orders across the mesh cells are also proposed. Then, to further extend this method, conforming tetrahedral meshes are replaced with non-conforming hybrid meshes. A modified DGTD formulation on hybrid tetrahedral/hexahedral meshes with p-adaptivity is presented in section 3.2. Then, we present the main details about the generation of hybrid meshes, especially for PV structures. Finally, in section 3.3, we show the impact of the mesh hybridization with p-adaptivity on the convergence of the scheme and on the computational cost of the simulations.

Locally adaptive DGTD method

In this section, the DG formulation presented in section 2.4 is modified to account for variable polynomial orders. Two distribution strategies are presented and compared. Resulting algorithms are tested for sequential executions, and the CPU and accuracy gains are evaluated on a realistic configuration relevant to photovoltaics.

Strategy based on the computed time-step per cell

We recall here the strategy initially introduced in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF], which is our starting and reference point. Starting with a given mesh, it seems evident that the final repartition of interpolation orders across the different cells will have a significant impact on the obtained accuracy, as well as on the CPU time required to obtain the numerical solution. Suppose that the solution is obtained on the given mesh with a homogeneous polynomial order P k . The point is here to see how, with a proper distribution of polynomial orders P l with l ≤ k, a solution with at least, a similar accuracy can be obtained at a lower computational effort. We start by defining the mesh size heterogeneity quantity as follow:

h = h max h min .
We also need to define the time-step in each cell T i as ∆t i , which is computed following the formula given in section 2.5.3, i.e.

∆t i = CFL(p)h i .
while ∆t p i represents the effective time-step obtained if cell T i is discretized with a P p polynomial expansion:

∆t p i = CFL (p)∆t i ,
with the second CFL factor CFL equal to:

CFL (p) = CFL(p) CFL(k)
with k the initial interpolation order. We also define the normalized time-step, which includes the computational load induced by the polynomial order, roughly estimated by the number of degrees of freedom (d.o.f.):

∆t p i = ∆t p i n(p) = CF L (p)∆t i n(p) ,
where n(p) is the number of d.o.f. in a P p cell. Finally, we define p min and p max the minimal and maximal user-authorized orders. In comparison with what was done in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF], we do not restrict the order jump between two adjacent cells (in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF], a P k -P l interface needs to fulfill the condition |k -l| = 1), even if it means that we should store more matrices in memory.

The first step of the repartition procedure consists of computing all the local ∆t i , and sorting them by ascending order. The cell with a lower time-step receives the minimal order p min , and we compute its normalized time-step ∆t p min 1 . We define two temporary variables, p loc and ∆t loc , which are respectively the current order assigned to elements and the current normalized timestep. For a given cell, the normalized time step for increased order p loc + 1 is compared to the current limiting normalized time-step ∆t loc . In the case where the first is higher than the second, switching to higher order is assumed to have a limited impact on the final time-step. Hence, p loc is increased by one, and ∆t loc is updated. The procedure is summarized in the algorithm 1.

Algorithm 1 Polynomial order distribution with the time-step strategy

1: for i ← 1, n t do
Compute time-step for each cell 2:

Compute ∆t i 3: end for end if 21: end for

Strategy based on the number of points per wavelength

It is well known that a proper numerical approximation of waves by a numerical method is related to the number of points per wavelength. It is generally accepted that 6 to 10 points per wavelength are usually sufficient to reach an acceptable accuracy.

In the context of this thesis, the main application is PV and, in general, simulations consist of the illumination of a device by a polychromatic plane wave over the whole visible spectrum [300 nm, 1000 nm]. In addition to this broad spectrum, the heterogeneity of solar cells, coupled with the complexity of the material used plays an essential role in the light propagation. Of course, when light travels from one medium to another, its speed changes, and so does its wavelength. In order to obtain accurate results and capture all the resonances, some physical conditions have to be fulfilled.

We start by defining in each cell T i , the metric d dof (i), which represents the maximum distance between two consecutive degrees of freedom. Then for the particular material in T i , the minimum effective wavelength λ(i) is computed. This is done by computing the wave's speed inside the material for all the wavelengths in the spectrum of interest. From there, it is possible to compute the minimum number of points per wavelength inside the tetrahedron T i as:

n ppw (i) = λ(i) d dof (i) . (3.1)
Then, given p min and p max the minimal and maximal user-authorized orders, and n ppw the minimum number of points per wavelength required, we are able to distribute the polynomial orders accross T h . We start by setting p min in all cell T i . Then we compute d dof (i) and n ppw (i). If n ppw (i) ≤ n ppw (i) and p(i) < p max , we increase p i . The procedure is summarized in algorithm 2.

Algorithm 2 Polynomial order distribution with the Points per wavelength strategy 1: p(:) ← p min Each cell receives the minimal order 2: for i ← 1, n t do Compute maximal distance between two d.o.f. for each cell 3:

Compute d dof (i) 4:
Compute n ppw (i)

5: if n ppw (i) > n ppw then 6:
Exit One can take a coarser mesh in this area

7:
end if

8: if n ppw (i) < n ppw then 9: if p(i) + 1 ≤ p max then 10: p(i) ← p(i) + 1 11: Cycle 12:
else 13:

Exit

One should refine the mesh here or increase p max 14:

end if 15:
end if 16: end for In the next part, we show on one example, the behavior of each strategy, and in particular that the second one performs better in the context of numerical simulation of PV cells.

Performance assessment

The first algorithm has been assessed on several nanophotonic cases in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF]. In this work, we want to highlight the main differences between the two strategies, and especially the configurations in which the second one is more suitable. As a reminder, it is shown in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF] that for a similar accuracy, a computational speedup of 6.15 for a locally refined cubic cavity, 2.60 for a metallic nanolens, and 2.19 for a bowtie nanoantenna is achieved when using local order polynomials. The second strategy is more based on physics than on numerics. Of course, in the end, the objective is still the same, i.e. making the numerical simulation cheaper, both in terms of memory and of computational time. However, this should be done such that the obtained results are accurate enough. For example, we consider a nanostructured PV cell composed of two gratings [WYL + 12]. This structure is studied in details in section 4.5.1 of chapter 4. The structured layer is shown in figure 3.1 with colors representing a metric based on the volume of each cell (blue cells are the largest ones; red cells are the smallest ones).

As one can see, some areas are critical, mostly at the top and bottom of the cones. The repartition of the 34544 cells is shown in figure 3.2(b), for the strategy based on the time-step and for the strategy based on the number of points per wavelength with n ppw = 6. We see that the first strategy in 3.2(a) is mostly leading to a P 2 interpolation everywhere, with a few elements using P 1 and P 3 interpolations. Regarding the second strategy, the repartition is shown in 3.2(b) for n ppw = 6. We note that the repartition is not exactly as shown, i.e. proportions are correct, but elements with P 1 interpolation are not only the smallest ones. Of course, the interpolation order assigned to a cell depends on the size of the cell, but also on the material and on the considered physical problem. We observe that most of the elements are rather using P 1 and P 4 interpolations. In figure 3.3, we compare both strategies on the computation of one important quantity of interest in the PV field: the reflection (noted R and defined in chapter 4). Both strategies are plotted against a converged DGTD-P 4 simulation. In figure 3.3(a), we observe a notable difference between DGTD-P 4 and DGTD-P 1 -P 4 (actually DGTD-P 1 -P 3 ), especially on the first part of the spectrum [200 nm, 320 nm], which corresponds to the shortest wavelength. The results of the second strategy are shown in figure 3.3(b). Three different values of n ppw are tested: n ppw = {4, 5, 6}. We observe a better behavior than the first strategy, even for 4 points per wavelength. Finally, computation times are given in table 3.1. As expected, the time-step based strategy outperforms the other ones in terms of computational time, with a speedup of 4. However, the solution is not good. For the second strategy, the speedup is directly related to the number of points per wavelength required. We obtained a speedup which goes from 3 for n ppw = 4 to 1.95 for n ppw = 6. To conclude this section, the strategy to be adopted depends mostly on the accuracy expected. When speedup is more important than accuracy, the first strategy and the second strategy with a small number of points per wavelength are preferred. When accuracy is still the main expectation, the second method with n ppw ≥ 5 should be preferred.

DGTD on hybrid meshes

As seen in the previous section, in order to be able to reach a certain accuracy, one should refine the mesh or increase the polynomial interpolation order in suitable regions of the computational domain. Each of these possibilities drastically increases the computing resources required for the simulation. However, it is common knowledge that hexahedral meshes are economical in terms of space (one hexahedron roughly replaces six tetrahedra) and advantageous in terms of accuracy. Here, we adapt the strategy presented in [START_REF] Durochat | Méthode de type galerkin discontinu en maillages multi-éléments pour la résolution numérique des équations de maxwell instationnaires[END_REF], which consists of combining a structured hexahedral mesh for the light propagation in untextured and homogeneous regions (here thick layers) with a fully unstructured tetrahedral mesh for the propagation in the nanotextured interfaces. Such meshes can be globally conforming or non-conforming and are naturally handled by a DG formulation.

Hybrid configuration

For continuous Finite Element Method (FEM), it is mandatory to use conforming meshes because of the continuity of the functions on the interfaces. Conforming hybrid meshes in 3D can be used by considering so-called transition elements like prisms or pyramids, which is not obvious from a meshing point of view. However DG methods, thanks to the absence of continuity of basis functions at the cells interfaces, can naturally handle hybrid and non-conforming meshes without a full modification of the core of the formulation presented in chapter 2.

For the sake of simplicity, here, the studied configurations are limited to a particular type of non-conformity. The use of non-conformal and hybrid faces is restrained to the case where neighboring hybrid faces match precisely, as shown in figure 3.4. In addition, the polynomial interpolation order in the parent cell have to be higher than the polynomial interpolation order of the small cells (see figure 3.5).

DGTD on hybrid meshes

To ease the reading, we consider the non-dispersive time-domain Maxwell equations in this section. However, the extension of the formulation is presented in equation (3.8) to the case of Maxwell's equations for dispersive media. Let Ω ⊂ R 3 be a bounded convex domain, and n the unitary outward normal to its boundary ∂Ω. The continuous, normalized Maxwell problem is: find

(E, H) in (L 2 (Ω)) 3 × (L 2 (Ω)) 3 ) µ r ∂H ∂t = -∇ × E, ε r ∂E ∂t = ∇ × H -J,
We consider Ω h as a discretization of Ω, relying on two quasi-uniform tesselations T h and

H h verifying T h = N T i=1 T i and H h = N H i=1 H i , where N ξ ∈ N *
is the number of mesh elements in each part with ξ ∈ {H, T }, and (T i ) i∈ 1,N T and (H i ) i∈ 1,N H the set of simplices. For i ∈ 1, N ξ , let C i be a tetrahedron or a hexahedron of the hybrid mesh. The internal faces of the discretization are denoted a ik = C i C k if C i and C k are neighboring cells, and n ik is defined as the unit normal vector to the face a ik , oriented from

C i toward C k . For each cell C i , V i is the set of indices {k ∈ 1, N | C i C k is a face}.
Let us define V h as:

V h = v ∈ L 2 (Ω) 3 , v T i ∈ (P p (T i )) 3 ∀ T i ∈ T h v H i ∈ (Q q (H i )) 3 ∀ H i ∈ H h
where P p (T i ) is the space of polynomials of maximum degree p on T i and Q q (H i ) is the space of polynomials of maximum degree q on H j . The semi-discrete fields, sought in space V h , are hereafter denoted (H h , E h , J h ), and on each cell T i the restrictions

(H i , E i , J i ) = H h T i , E h T i , J h T i are defined. A set of scalar basis functions (φ ik ) 1 k d i is defined for each T i , where d T i is the number of degrees of freedom (d.o.f.) per dimension (respectively (ψ jk ) 1 k d H i
for each H i ). Additionally, to each scalar basis function, the three vectors φ v ik are associated (respectively ψ v ik ):

φ x ik =   φ ik 0 0   , φ y ik =   0 φ ik 0   , φ z ik =   0 0 φ ik   , ψ x ik =   ψ ik 0 0   , ψ y ik =   0 ψ ik 0   , ψ z ik =   0 0 ψ ik   ,
One now seeks the approximations E h of E in space V h (respectively H h of H). The contribution of each cell is therefore defined as E i = E h T i . Here, one must notice that, for a 3D system, E i is actually a vector that has 3 components:

E i =   E x i E y i E z i   ,
each of which is locally expanded on the chosen set of basis functions:

E v i = d i j=1 E v ij φ ij , v ∈ {x, y, z} . (3.2)
Therefore, for practical purpose, one defines three vectors of d i components:

E v i =    E v i1 . . . E v id i    , v ∈ {x, y, z} ,
as well as the following 3d i components vector:

E i =       E x ij 1 j d i E y ij 1 j d i E z ij 1 j d i       .
The same four steps are done for H and for an element H i , even if it is not detailed here. We consider the weak formulation of Maxwell's equations without current:

ˆCi µ r ∂H ∂t • φ v ik + ˆCi E • ∇ × φ v ik = l∈V i ˆail (E * × n il ) • φ v ik , ˆCi ε r ∂E ∂t • φ v ik - ˆCi H • ∇ × φ v ik = - l∈V i ˆail (H * × n il ) • φ v ik . (3.3)
We can develop the first volumic integral using the local field expansions (3.2):

ˆCi ε r ∂E i ∂t • φ x ik = M εr i ∂ E x i ∂t k
as well as the curl integral:

ˆCi H i • ∇ × φ x ik = -K i × H i x k ,
without any changes with respect to the formulation presented in chapter 2. Indeed, the main changes will appear in the evaluation of the flux integral on hybrid faces. For the sake of simplicity, we consider the centered flux. The generalization to the upwind flux is straightforward. Consider a neighbor cell H l of T i with a Q q polynomial expansion. The xcomponent of the flux integral on their common face a il from the T i side is:

ˆail (H * × n il ) • φ x ik = ˆail H y * n z il -H z * n y il φ ik = ˆail H y i + H y l 2 n z il - H z i + H z l 2 n y il φ ik (3.4)
Consider the following expansions for H y i and H y l , and the analogous relations for H z i and H z l :

H y i = p q=1 H y iq φ iq and H y l = q r=1
H y lr ψ lr .

(3.5) Plugging (3.5) in (3.4) leads to:

ˆail (H * × n il ) • φ x ik = ˆail n z il 2   p q=1 H y iq φ iq φ ik + q r=1 H y lr ψ lr φ ik   - ˆail n y il 2   p q=1 H z iq φ iq φ ik + q r=1 H z lr ψ lr φ ik   , = 1 2 p q=1 H y iq n z il -H z iq n y il ˆail φ iq φ ik + 1 2 q r=1 H y lr n z il -H z lr n y il ˆail ψ lr φ ik , = 1 2 p q=1 H iq × n il ˆail φ iq φ ik + 1 2 q r=1 H lr × n il ˆail ψ lr φ ik , = S il H * ,i × n il x k + S * il . H * ,l × n il x k , (3.6)
where, in accordance with the definitions of chapter 2, the flux matrices are:

(S il ) jk = ˆail φ ij φ ik and (S * il ) rk = ˆail ψ lr φ ik . (3.7)
To summarize, the flux integrals are cut in two parts: (i) the part corresponding to local contribution, which is integrated via the "regular" flux matrix, and (ii) the part corresponding to the neighbor contribution, which is integrated via non-conforming hybrid matrices. For a P p -Q q interface, the flux matrix is rectangular, of size s(p) × s(q) (we recall that s(p) = (p+1)(p+2) 2 for a tetrahedron and s(q) = (q + 1) 2 for a hexahedron). From this point, the rest of the derivation of the numerical form of the scheme is similar to the tetrahedral mesh case, and the reader is referred to section 2.4 for details. The final semi-discrete scheme is:

M µr i ∂ H i ∂t = -K i × E i + l∈V i S il E * ,i × n il + l∈V i S * il E * ,l × n il , M εr i ∂ E i ∂t = K i × H i - l∈V i S il H * ,i × n il - l∈V i S * il H * ,l × n il , (3.8)
where the definition of S * il is easily deduced from (3.7) and the analogous definition of S il given in section 2.4. Time integration still makes use of the schemes presented in section 2.5.

Numerical study of h-convergence

In this section, we study the numerical convergence of the DGTD method on hybrid meshes. For this purpose, we consider the cubic cavity case with PEC boundaries presented in chapter 2. Four hybrid meshes with increasing refinement are used. The cubic domain is split into two parts, as shown in figure 3.6, one for each element type. Upwind fluxes and the optimized 14 stages fourthorder low-storage Runge-Kutta time scheme developed in [START_REF] Niegemann | Efficient low-storage runge-kutta schemes with optimized stability regions[END_REF] are used. The obtained results are given in table 3.2. As expected for a P p -Q q configuration, the h-convergence rate (r) is limited by the lowest approximation order which is min(p, q) + 1.

P p Q q (a) Test configuration for h- convergence (b) M4

Numerical study of hybrid interfaces

In a second step, we want to highlight the potential presence of parasitic reflection near the hybrid interfaces. For this purpose, we consider a cubic domain with absorbing boundary conditions, with the injection of a polychromatic incident plane wave in the x direction. This choice is motivated by the simulations of light trapping in PV cells, which in most cases uses this type of source. For the mesh, we construct a sample of possible non-conformity (see figure 3.7). All meshes are based on a cartesian grid of 216 hexahedra (6 × 6 × 6). Tetrahedral zones are obtained by splitting hexahedra in 6 tetrahedra. In figure 3.8(a), we compare the error on a tetrahedral mesh (DGTD-P 1 ), with the two examples from the first configuration: DGTD-P p -Q q means we inject the plane wave in the tetrahedral region (P first), while DGTD-Q q -P p means we inject the plane wave in the hexahedral region (Q first). In figure 3.8(b), we consider the case where the hybrid interface is splitting the injection surface in half. For both graphs, the error is similar, with a slight increase for the second configuration.

M1 M2 M3 M4 E -E h r E -E h r E -E h r E -E h r P1 -Q1 5.05 × 10 -2 -1.

Realistic PV cells

PV structures present a large variety of configurations. However, most of them share common properties. Indeed, they are generally constituted of several stacked layers of materials, and the interface between each of them can be nanostructured for light trapping purposes [Pae12]-[WYL + 12]. However, those nanotextures are only present at the junction between two layers and represent a very thin portion of the solar cell (a few nanometers compared to several hundreds, as shown in figure 3.9). Therefore, the use of hybrid meshes is relevant in this context. The generation of hybrid meshes is possible by using the Gmsh meshing software [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] 1 . Even if it is not really a proper feature from Gmsh, its ease and modularity make it possible to generate non-compliant and hybrid meshes of PV stacks quite easily. The main idea is to build the hybrid interfaces explicitly in the Gmsh script. This means that we must select in advance the regions that are meshed with hexahedra. This corresponds to having a set of intervals [h i min , h i max ], i ∈ N * , which represents the delimitation of each area.

The main restriction with this approach is that each hexahedral zone must have the form of a block. Our second possibility consists in using an intermediate tool, developed initially to study the radar cross section (RCS) of aircraft. This tool allows us to take as input a surface or volume, tetrahedral or hexahedral (or already hybrid) mesh, as well as a domain size (l x × l y × l z ), a max hexahedron size and a minimum hexahedron size. From these information, the software builds the hybrid interface and the hexahedral zones around it. If the input meshes are surface, the Gmsh tool is used to fill the gap between interface and surface.

P p Q q (a) Config. 1 Q q P p (b) Config. 2 Q q P p P p Q q (c) Config. 3 Q q P p (d) Config. 4

Numerical assessment for a double conical grating cell

We consider the case of the nanostructured PV cell with the two conical gratings again [WYL + 12]. The homogeneous regions are removed and replaced with layers of hexahedra. It concerns the silicon layer surrounded by the gratings, but also the artificial ones (PMLs and SF areas). One can see the hybrid mesh and its corresponding element's size in figure 3.10. For this setup, the same pattern as in section 3.1.3 is observed: the largest elements are hexahedra placed between the two gratings, while the finest ones are tetrahedra located at the extremity of cones.

Results are shown in figure 3.12. We start from the best working configuration using the point per wavelength strategy, as shown in part 3.1.3, i.e. using n ppw = 6. The results are matching on the range [230 nm, 500 nm] for both P 1 -P 4 -Q 3 and P 1 -P 4 -Q 4 interpolation orders. However, for the very first wavelengths, there is a lack of reflection when Q 3 interpolation is used. We also plot in figure 3.11 the contour lines of the amplitude of the DFT of E at λ = 450 nm.

The computational times are shown in table 3.3. The results for the original DGTD-P 4 method and the improved DGTD-P 1 -P 4 version using local interpolation order with n ppw = 6 are recalled, as well as the two hybrid configurations considered. The speedup in parentheses corresponds to the speedup from the DGTD-P 1 -P 4 simulation. The DGTD-P 1 -P 4 -Q 4 , n ppw = 6 presents a speedup of 1.40 as compared to the p-adaptive simulation. We do not see a real overhead due to the construction and storage of the specific hybrid matrices. In any case, this overhead is mainly present in the initialization phase, which is of less importance compared to the time loop.

Numerical assessment for multilayer PV type nanostructures

Here we consider an idealistic PV stack composed of six layers. The interface geometry is manually defined using a(x, y) = cos y 2 , which gives a bi-periodic domain. We consider five different For this problem, the quantity of interest considered is the volumetric absorption (defined in chapter 4). A convergence analysis based on this quantity is done in figure 3.14. The volumetric absorption profile obtained with the coarse hexahedral mesh has a general appearance close to the reference solution. However, resonances near 500 nm are not well captured. On the contrary, the primary resonance is not in the same place. It is important to note that this anomaly does not come from a convergence concern but from a geometric representation error due to the staircasing effect. The refinement of the cartesian grid has been done in such a way as to obtain an equivalent solution. It results in a hexahedral mesh of 320000 elements. Mostly because of the simplicity of the textured interfaces, we obtained almost the same results for both tetrahedral meshes. The hybrid mesh generated is using two levels of hexahedron refinement in order to restrict the tetrahedra area. Due to the high ratio between the smaller and the lower, we obtained better results using a P 2 interpolation for the tetrahedral region while using a mixed Q 2 -Q 3 for the hexahedral one. Subfigure 3.14(f) show the comparison of each converged results, which are all really close, except for the coarse hexahedral mesh. Now that we see that M2, M3, M4 and M5 can achieve the same (good) result, we want to highlight the computation cost of all of them. The simulations have been performed sequentially on an Intel(R) E5-465L (2.60 GHz -8 cores) with 512 go of RAM. We choose this workstation due to its huge memory, which allows us to compute our quantity of interest on M2 and M4. The computation times are shown in table 3.5.

A tot (λ) DGTD-P1 -Q1 DGTD-P1 -Q2 DGTD-P1 -Q2 -Q3 DGTD-P2 -Q2 -Q3 (e) M5
A tot (λ) M1 M2 M3 M4 M5 (f) Comparison

Conclusion

In this part, we have presented two numerical strategies that enable a better resolution of complex geometrical problems. In the first part, we recalled some preliminary work done by Viquerat et al in [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF] about local polynomial approximation. In this work, we extended the latter with a criterion based on the number of degrees of freedom per wavelength. In comparison with [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF], this approach allows us to have an a posteriori control of the behavior of the simulation from a starting mesh. Some interesting results have been shown, with a number of points per wavelength between 6 and 8. This allows an easier use of very various mesh sizes. In a second time, we presented the modifications required in the DGTD formulation in order to take into account hybrid meshes.

A convergence analysis has been performed on a basic test case, and we investigated the behavior of several hybrid interfaces. From our examples, we did not experiment with any notable degradation with an incident normal polychromatic plane wave. Finally, in the last part, we introduced a multilayer PV structure. We highlighted the influence of the staircasing effect, which is one of the major drawbacks of the use of cartesian grids. Then, performances of the previous methods were compared, with a clear improvement when combining hybrid non conform meshes with local polynomial interpolation, both in computational time and in number of d.o.f.. Photovoltaics (PVs) is a very prominent field of research, with widespread awareness of renewable energy. One of the leading scientific concerns around PVs lies in the discovery of processes and strategies to improve the performance of solar cells, in order to enable them to become the reference green energy. Since the beginning of PVs, silicon has been the primary material used as an absorber, despite its significant defects. This is majorly caused by its low price. One of the main way to overcome these failures is to use light trapping to focus the light in the absorber. Light trapping consists in structuring one (or more) layers so that the light is trapped in the absorbing material. This structuring is the subject of much research [WYL + 12], [START_REF] Eyderman | Near perfect solar absorption in ultra-thin-film {GaAs} phtonic crystals[END_REF], [START_REF] Lockau | Optical modeling of thin film silicon solar cells with random and periodic light management textures[END_REF]. Knowing that the costs and time required to manufacture the cell are not negligible, the use of numerical simulations seems to be an asset and a primary possibility for the testing and calibration of sophisticated PV devices.

The main criterion to evaluate and sort solar cells is the efficiency. It corresponds to the percentage of light received and transformed into electrical power. This transformation involves several physical processes, including optics and electricity. In the context of this thesis, we only consider the optical part, which consists in quantifying the proportion of incident solar light absorbed in a specific material, generally a semiconductor. We will not deal with the generation of electron hole pairs, recombination, and carrier transport, which are the other fundamental parts of solar cell.

Three possibilities occur when light reaches the solar cell: (i) the incident photons are reflected at the top surface, (ii) they are absorbed in the material, or (iii) the light is transmitted through the cell. For most of the PV devices, both reflection and transmission are considered as loss given that photons are not absorbed, and therefore, do not generate power. Both issues can be taken into account in different ways. Usually, reflection is handled using anti reflective coating [START_REF] Raut | Anti-reflective coatings: A critical[END_REF]. This technology is not restrained to PV; it is also used in a wide variety of applications where low reflection is desired, e.g. on corrective lenses, and camera lens elements. To overcome the transmission, a mirror is often used at the back of the cell, reflecting the light initially transmitted and allowing a second pass through the absorber. This allows to exceed the single pass limit absorption. In recent research, those two strategies have been used in order to improve PV cell efficiencies. For example, recently in [CCDL + 19], Chen et al used a nanostructured back mirror to reach the efficiency record for ultrathin solar cells (nearly 20%). The simulation of this part can be done by solving full-wave Maxwell's equations.

In this chapter, we first review the most commonly used numerical approaches for the simulation of light trapping in solar cell devices in section 4.1. In section 4.2, we recall the definition of the figures of merit that we will use for assessing the efficiency of such devices from the point of view of solar light absorption, such as the reflection and the transmission coefficients, External Quantum Efficiency (EQE) and the short circuit current density (Jsc). Then we present the similarities shared by PV geometries that we usually encounter, namely their bi-periodicity. We also present some methods to manage high definition measured textures data with noise. Those challenging data are preprocessed to create a good and usable mesh for our DGTD solver. We discuss both meshing aspects, as well as the observable of interest, specific to PV applications in section 4.3.

All these elements are illustrated with several realistic test cases from PV literature or studies undertaken in collaboration with physicists in the PV field in section 4.4, 4.5 and 4.6.

Numerical modeling of light trapping in PV devices

Solar cell physics is composed of several branches. The two main ones are: (i) optics, which characterize the interactions between solar light and the solar cell and, (ii) semiconductor physics, which deals with the generation of electric current from carrier generation and transport. In order to simulate each of those parts, it is necessary to use different methods, each adapted to model the physics of interest. For the semiconductor physics, one can solve the drift-diffusion model [ACM].

For the optical part, which is the physics we are considering in this thesis, several methods exist for the computation of solar cell properties. To cite a few of them, the Ray Tracing method is well adapted for textures with characteristic sizes that are larger than the relevant wavelength [START_REF] Baker-Finch | Reflection of normally incident light from silicon solar cells with pyramidal texture[END_REF]. This method literally consists in tracing a ray of light using simple geometrical optics based on the refractive index of each material. If layers are thin but still unstructured, it is possible to use the Transfert Matrix (TM) method. TM method is based on continuity conditions for the electric field between two media. Knowing the field at one side of a layer allows us to compute the field at the other side from a simple matrix operation. The Rigorous Coupled-Wave Analysis [START_REF] Moharam | Rigorous coupled-wave analysis of planar-grating diffraction[END_REF] is suitable for solving scattering from relatively small periodic structures. When texturization details are smaller or close to the relevant wavelength, it is important to take into account the optical wave effects. This is possible via discretization methods, that solve electromagnetic waves equations such as the FDTD [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF] or the Finite Element Frequency-Domain (FEM) [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF] methods. Recently, some coupling methods have emerged. They proposed to couple different simulators (usually RCWA, ray tracing, and TM) in order to combine the benefits of each. The two main methods are OPTOS [TEK + 15] and the Coupled Modelling Approach (CMA) [TSL + 15].

Rigorous Coupled-Wave Analysis (RCWA)

The Rigorous Coupled-Wave Analysis, introduced in 1981 by Moharam [START_REF] Moharam | Rigorous coupled-wave analysis of planar-grating diffraction[END_REF], is a frequencybased, semi-analytical optical simulation method. It is mostly used to solve scattering from periodic dielectric structures. The RCWA method starts by discretizing the domain into layers in the z direction. Curved gratings need stairstep-approximated layers (see figure 4.1). Electromagnetic modes are calculated in each layer and then analytically propagated through the layers by matching boundary conditions at layer interfaces, in the same fashion as TM. These modes are expanded using Fourier series, which leads to an infinitely large system of Ordinary Diffential Equations (ODE). A finite-dimensional problem can be reached by truncating the Fourier expansions, at the cost of the accuracy of the method. RCWA method is very efficient for the calculation of far-field reflection and transmission coefficients, but the accurate computation of the near field at material boundaries remains challenging. RCWA has been used for many optical simulation of nanotextured solar cells [ABK + 19],[BEH + 16], [ACM]. The main drawback of the RCWA method is that it is not suited for non-cartesian geometries (see figure 4.1). Moreover, even if the method is dedicated to periodic structures, the number of modes necessary can become enormous. Indeed, it is related to the size of the elementary pattern and of the relevant wavelength. This generally leads to huge matrices to be inverted.

Finite Difference Time-domain (FDTD)

The Finite Difference Time-domain (FDTD) method is undoubtedly one of the most popular methods for time-domain computational nanophotonics. This method is also widely used for the study of light trapping in solar cell devices. The method relies on two Taylor expansions for both spatial and temporal derivatives. The popular version presented by K.S. Yee in 1966 [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF] relies on a particular discretization of the computational space. The spatial element used is now referred to as Yee cells (see figure 4.2(a)). In the original Yee scheme, the spatial derivatives are discretized using second-order central differences while time integration is achieved with a secondorder leap-frog scheme. As of today, FDTD represents an easy to use and efficient method to solve electromagnetics problems, combining simple implementation and high computational efficiency. However, the FDTD method also presents some limitations. First, a smooth discretization of curved geometries is impossible due to the fixed cartesian grid imposed by the Yee algorithm. This approximation leads to the well-known staircasing effect (see figure 4.2(b)), which is an important source of inaccuracy [DDH]. To overcome this pitfall, one can either use an extreme refinement of the grid, which leads to a drastic rise of the computational cost, or exploit one of the numerous possible modifications of the FD method that have been proposed for tackling the staircasing effect [START_REF] Hao | Analyzing electromagnetic structures with curved boundaries on cartesian (fdtd) meshes[END_REF]. However, these modifications represent a tradeoff between the simplicity of the classical algorithm and the accuracy of the boundary description. An other source of inaccuracy in the FDTD method is faced when dealing in the case of heterogeneous problems. Indeed, since the electromagnetic fields are not smooth at the interface of different medium, the Taylor approximations used are no longer valid. More advanced FDTD methods aim to address this issue [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]. Unfortunately this is once again at the price of an increased complexity of the algorithm.

Finite Element Method (FEM)

Continuous FEM is used as an alternative method to the FDTD method. The FEM for CEM was introduced in 1969 by Silvester to solve waveguide problems [START_REF] Silvester | Finite element solution of homogeneous waveguide problems[END_REF]. Unlike the FDTD, here the domain tessellation is an unstructured tetrahedral or hexahedral mesh. A discrete variational formulation is obtained from the continuous Maxwell's equations by approximating the unknowns in a finite-dimensional space. Then, the discretization of this weak form leads to a sparse system of equations in terms of the whole set of d.o.f.. This system is one of the drawbacks of FEM methods in time-domain since it must be solved at each time-step. Moreover, one should avoid the use of nodal basis functions as it was proved that they could lead to spurious modes, due to an ill representation of the curl kernel [START_REF] Sun | Spurious modes in finite element methods[END_REF]. Edge finite elements (also known as Nédélec edge elements) were introduced in 1986 by Nédélec [N 80]. These elements display several interesting properties: (i) their divergence is zero, and (ii) each basis function associated with an edge has a constant tangential component on the latter, and a zero tangential component on the others (see figure 4.3). This leads to a naturally enforced tangential continuity of the electric field across the edges. In order to adjust the accuracy of the simulation, FE methods can use either (i) a local refinement or coarsening of the mesh, (ii) a local or global increase of the order of the basis functions (if polynomial, it corresponds to an higher polynomial order), or (iii) a combination of both. However, these improvements lead to larger linear systems to solve at each time-step, which can make the FE method impractical in time-domain simulations for huge systems. For this reason, FE methods are more often used in frequency-domain. However, a few references can be found exploiting time-domain FE for nanophotonics applications [START_REF] Huang | Modeling backward wave propagation in metamaterials by the finite element time-domain method[END_REF].

Physical quantities of interest

When assessing the quality of a PV device, some fundamental quantities allow to measure the performances of different criteria such as its anti-reflection capability or the maximum current than can be generated from a solar cell. In this section, we introduce the figures of merit that we will use in our simulations.

Reflection, transmission and absorption

First of all, we need to introduce the Poynting vector. It is an appropriate quantity to describe the energy flux of a propagating electromagnetic wave. Its definition in time-domain is:

π (x, t) = E (x, t) × H (x, t) . (4.1)
More precisely, we are interested in the time-averaged Poynting vector which is define as:

π (x, λ) = 1 2 E (x, λ) × H * (x, λ) , (4.2)
with E and H the Fourier transform of respectively E and H and λ the wavelength. Given that definition, consider the situation of figure 4.4, where the scatterer is enclosed between two TF/SF surfaces. Between the two surfaces, i.e. in the total field region, the computed fields are E tot and H tot and one as:

E tot = E inc + E sca , H tot = H inc + H sca .
The incident field (E inc , H inc ) is imposed on the TF/SF interface, and the computed fields in the scattered region are E sca and H sca . Using these fields we can define π inc and π sca as:

π inc (x, λ) = 1 2 E inc (x, λ) × H * inc (x, λ) , (4.3)
π sca (x, λ) = 1 2 E sca (x, λ) × H * sca (x, λ) , (4.4)
which are the key ingredients for the calculation of the physical quantities of interest. Now we can define the reflection (R) and transmission (T ) coefficients, which are two wavelength-dependent observables. When a scatterer is illuminated by an incident field, R represents the proportion of energy reflected to the injection plan, while T is the amount of energy that traveled through the system: and:

R (λ) = ˆSin π sca (•, λ) • n ˆSin π inc (•, λ) • n , ( 4
T (λ) = ˆSout π tot (•, λ) • n ˆSin π inc (•, λ) • n , (4.6)
where S in and S out are the injection and output planes defined in figure 4.4 and n is the outgoing normal vector. Additionally, we can define the absorption (A) from R and T as:

A (λ) = 1 -T (λ) -R (λ) . (4.7)
which represents the amount of energy absorbed by the device. It is one of the most important quantity in the PV field as it represents the amount of energy absorbed by the cell, and so potentially transformed to electricity. The equation (4.7) point out two primary strategies to optimize solar cells absorption: (i) one can use an anti-reflection (AR) coating at the front surface of the cell to reduce R, (ii) or one can use a backside mirror, most of the time made of silver, to minimize T and lead to light trapping.

Remark : For non-periodic scattering problems, the absorption cross sections (C abs ) have to be used in order to compute the absorption [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF].

In the following, A computed from (4.7) will be noted A tot since it corresponds to the total absorption in the PV cell.

Volumetric absorption

Absorption can also be computed with a volumetric method. It is possible to evaluate the Ohmic losses (W) inside a given material instead of computing the flux of the Poynting vector through surfaces. It is shown in [START_REF] Landau | Electrodynamics of Continuous Media[END_REF] that the power absorbed by the scatterer as Ohmic losses is:

A s (λ) = P Ohm (λ) = ε 0 λ 2 ˆΩs (ε r (λ)) E(r, λ) 2 , (4.8)
where Ω s is the volume of the scatterer s. If the domain is composed of only one absorber, one can choose which quantity to compute, i.e. the volumic or the surfacic absorption. As discussed in [START_REF] Andreas | Pitfalls in the Determination of Optical Cross Sections From Surface Integral Equation Simulations[END_REF], the computation of the absorption via the Poynting vector integration shows large errors when there is a strong scattering and can even produce unphysical negative absorption. These effects are well explained in the paper mentioned above. However, evaluation of the volumetric absorption is way more expensive as it requires to compute an integral over a volume rather than a surface.

Quantum Efficiency

Quantum Efficiency (QE) or Incident Photon to Converted Electron (IPCE) is a frequency-dependent quantity used to measure the efficiency of a PV device. It is one of the key quantities as it gives information on the current that a given cell will produce when illuminated by a particular wavelength.

There are two quantum efficiency definitions used for the study of solar cell devices: the EQE, which is the ratio of the number of electrons generated to the number of incident photons per second and the Internal Quantum Efficiency (IQE), which is the ratio of the number of electrons generated to the number of incident photons absorbed by the cell per second as summarized in the following equations: Considering that the DGTD solver used in this thesis is only dealing with optical effects, we can assume that all absorbed photons are collected at the junction without recombination. Then, the EQE is computed as:

EQE = electrons/sec photons/sec , ( 4
EQE(λ) = P Ohm (λ) I 0 • S dom = A(λ), (4.11)
where I 0 is the intensity of the incident light and S dom is the surface area of the periodic pattern perpendicular to the incident light.

Remark : Unlike the absorption calculated from R and T , volumetric absorption can be calculated in various materials. However, the volumic absorption in each layer must coincide with A tot .

Short circuit current and short circuit current density

One other important quantity is the short circuit current (Isc) (mA). It represents the current through the solar cell when the solar cell is short-circuited. The short-circuit current is due to light-generated carriers being generated and collected. It is therefore the largest current from the solar cell that can be drawn. It is more common to use the Jsc (mA / cm 2 ) which is the current produced per unit cell area, rather than the Isc because it removes the dependence of the solar cell area. Jsc mostly depends on two ingredients: the incident light spectrum or spectral irradiance, standardized to the AM1.5G spectrum for most of terrestrial solar cell measurements (see figure 4.5); and the EQE. From those two ingredients we can compute it as:

Jsc = ˆqEQE(λ)AM1.5g(λ)dλ, (4.12)
with q is the elementary electronic charge (C). Remark : Some of the quantities defined above are sometimes referred to as near-field quantities and others as far-field quantities. It is important to notice that the numerical convergence of each field type requires different prerequisites, which will be the topic of a forthcoming part.

Dealing with random textures

In this section, we explain the different steps required to move from random textures given as AFM (Atomic Force Microscopy) data, to a geometrical model that can be used for the simulations with our DGTD solver. Such a geometrical model consists of a fully unstructured tetrahedral or hexahedral/tetrahedral mesh, which is processed using available mesh generation tools. The two tools used in this thesis are Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF], a free 3D finite element mesh generator with a built-in CAD engine, and MeshGems1 , a suite of meshing software components, integrated in most of the CAD vendors.

Texturized structures from AFM

The main tool used to inspect the shape of textured layers is Atomic Force Microscopy (AFM).

It is a very high resolution type microscopy capable of a resolution on the order of fractions of a nanometer. AFM outputs a discrete planar grid defined such that (x, y) ∈ [0, L] 2 with x = nδh, y = mδh, δh the resolution between two measured points and (n, m) ∈ [0, N ] 2 , N = L/δh. For all grid points, we also have a corresponding altitude noted a(x, y) > 0.

From there, we already have, in some ways, a quadrilateral surfacic mesh of the layer. However, most of the time, these data present abrupt jumps and spikes, which occur because of the way the AFM is used. In order to remove these defects, we choose to consider our set of points as an image, with the altitude a representing the grayscale of each point (x, y). Doing so, we can use image smoothing techniques coming from the digital image processing field to remove high-frequency content such as noise or edges. Two types of smoothing were tried here. First, we chose to use a box blur, also known as box linear filter. This method is easy to implement as it is basically a local average of each point. Mathematically, it consists of convolving the image with a normalized box filter:

a post (x, y) = K * a(x, y), (4.13)
with a post the filtered image and K the filter kernel used here. The kernel K represents the shape and the size of the neighborhood taken into account. In this work we choose to take a square symmetric kernel:

K = 1 n 2      1 1 • • • 1 1 1 1 . . . . . . . . . 1 • • • 1      n×n , (4.14)
with n > 0, n ∈ N. For example, we obtain the identity by taking n = 1. An illustration of the convolution for a 5 × 5 input dataset and a kernel with n = 3 is shown in figure 4.6. By specifying that we only consider kernels with odd n, we can define n = n-1 2 . Then, we can develop equation (4.13) as:

a post (x, y) = n k=-n n l=-n K(k + n , l + n )a(x + kδh, y + lδh) (4.15)
Remark : This implementation of discrete convolution requires some adjustment when the kernel moved to positions where it does not fit entirely within the image, e.g. at the edges and corners. It is possible to add some "ghost" pixels in order to normally apply the kernel.

The second choice is called the Gaussian smoothing. The changes occurs at the kernel level, where the values (1/N 2 ) are replaced with Gaussian values. The kernel is now represented as

K = (G(x i , y i )) (i,j)∈{1,...,N } 2 , (4.16)
with

G(x, y) = 1 2πσ 2 e -x 2 +y 2 2σ 2 , (4.17)
where σ represents the standard deviation. The advantage of this strategy, in comparison with the blur box, is the natural account of the distance between the source and the nearby values. Indeed, in this case, each pixel is still a weighted average of its neighboring values, but the weight is decreasing with the distance. The matrix 4.18 is an example of a 5 × 5 Gaussian kernel obtained using σ = 1.0. Obviously, the smoothing introduces discrepancies between the original model and the obtained surfacic mesh. It would be essential to check, as a future step, how we can control and minimize the impact of the smoothing concerning the standard deviation σ and the kernel size n. The first idea is to find a ratio between the resolution of the point cloud and the accuracy needed for the calculation.

K = 1 273       1 4
      . ( 4

Bi-periodic structures

Once the acquisition of surfacic meshes of each layer is made, we still lack the four lateral faces to have a closed domain. In particular, we have to impose Periodic Boundary Conditions (PBC) that allow to simulate artificially infinite mono-or bi-directional structures while considering only one elementary pattern. This is done by matching cells from a periodic boundary face with their neighbors on the opposite boundary of the domain. It is evident that the randomly textured nature of the geometries considered here does not naturally allow this matching to be applied.

In order to apply this principle to our realistic solar cell meshes, we explore two possibilities. The first one consists of building a mesh by gluing together four meshes obtained by symmetries. It is done in two steps, as shown in figure 4.9:

1. Considering a domain Ω = [x min , x max ] × [y min , y max ] × [z min , z max ], discretized in a mesh
containing the vertices v i = (x i , y i , z i ) with i ∈ {1, ..., N }. Let P + x be the plane defined by x = x max and p(v i ) be the orthogonal projection of v i on P +

x . From here, we can get a new mesh replacing the points with:

v i = 2p(v i ) -v i , i ∈ {1, ..., N }. (4.19)
It should be noted that elements keep the same connectivity, which allows a very fast execution of the operation. From here we can merge the two meshes, with the only difficulty being to not duplicate the nodes on P + x . 2. The second step is straightforward. It consists in repeating the same operations from step 1 taking into account all the vertices of T 2 h, discretization of

Ω 2 = [x min , 2x max -x min ] × [y min , y max ] × [z min , z max ],
with the plane P + y defined by y = y max . This method presents two main drawbacks. The first and biggest one is that the obtained mesh is four times larger than the starting one, both in terms of computational domain size than in terms of elements. Secondly, one of the layers may present a hole or a bump at the edge of the domain. Joining both parts at this point can create problematic geometrical singularities, which can drastically affect the simulation results.

Remark : One should have a look at the mesh to see if it is more interesting to use P +

x or P - x in (i), respectively P + y or P - y in (ii). The second option is more intrusive and needs to be done before the volumic mesh generation. It uses a so-called Tukey window function on each textured layer to level out each edge, as proposed in [START_REF] Jäger | Ambiguities in optical simulations of nanotextured thin-film solar cells using the finite-element method[END_REF]. The main idea behind this procedure comes from the evident simplicity of matching rectangular faces. To obtain such rectangular faces in x and y directions, we apply the Tukey window function w to all of our point:

a tk (x, y) = w(x) * w(y) * a post , (4.20)
In what follows, we use the function defined on [0, 1] by:

w(ξ) =              1 2 {1 + cos[π( 2ξ r -1)]} 0 ≤ ξ ≤ r 2 , 1 r 2 ≤ ξ ≤ 1 - r 2 1 2 {1 + cos[π( 2ξ r - 2 r -1)]} 1 - r 2 ≤ ξ ≤ 1,
with ξ ∈ {x, y} and r ∈ [0, 1] is a parameter managing how smooth the transition between the edges and the textures will be. One can see in figure 4.10, an analytical example which shows how the Tukey window works. We also show the obtained surface starting from a real dataset in figure 4.11.

It is crucial to notice that the value chosen for r impacts the quality of the obtained mesh. Taking r ≈ 0 will restrict the dumping area to a tiny area, creating sharp edges on the border, similar to the second drawback of the first method. On the other hand, having r ≈ 1 will cause a loss of details on a large area of the domain. In addition to that, the flat areas located at the boundary can lead to the creation of cavity modes that should usually not exist. Once this is done, we can build a volumic mesh using a mesh generation tool, taking care to specify the need for bi-periodicity.

GaAs benchmark

In order to assess the capabilities of the DGTD method introduced in chapter 2 for the simulation of light propagation in nanostructured solar cells, we participated in a benchmarking study proposed by Stéphane Collin at C2N. The main objective of this study is to compare the different numerical methods for the simulation of light absorption in nanostructured PV cells. Except for the DGTD method which is provided by us, the three other methods are the ones presented in section 4.1. For this purpose, a series of solar cell configurations based on Gallium Arsenide (GaAs) material of increased geometrical complexity has been specifically devised for this study. Solar cells considered are constituted of multiple slabs, and the silverback reflector is nanotextured with a non-symmetric pyramidal grating. The horizontal size of the elementary pattern is 600 × 600 nm 2 . The solar cell is illuminated from above and a semi-infinite slab of silver is used at the bottom. PMLs are used

Material models

The optical properties of the different materials that constitute the considered device structure have been fitted to the parameters of the generalized dispersion model [START_REF] Viquerat | Fitting experimental dispersion data with a simulated annealing method for nano-optics applications[END_REF] presented in chapter 2. Different numbers of first-order poles and second-order poles have been tested in order to obtain the best possible results. The fitted models obtained are ranked according to their error on the real part of the permittivity and on the imaginary part of the permittivity. From there, we choose a compromise between the number of poles and the accuracy of the model. Indeed, the addition of a pole is equivalent to adding an ODE (or two for a second-order pole) to the Maxwell-GDM equations. The obtained permittivity functions are plotted in figure 4.14. As can be seen, they are relatively well approximated in the wavelength range λ = [300, 1000] nm. Note that the SiO 2 is here modelled with a constant permittivity and without losses, i.e ε r = 2.1025 and ε i = 0.

Numerical convergence

Numerical convergence has been assessed for Structure 1 only, although in all cases we performed simulations with the DGTD method based on P 4 interpolation in addition to P 3 ones to ensure there were no major variations in the obtained results. As an example, we plot in figure 4.15 the total absorption obtained for Structure 1, as a function of the polynomial order. As can be seen, P 3 and P 4 results are very close, except for some small discrepancies in the high wavelength region of the spectrum. For the sake of completeness, we also compare the spectra obtained using the three meshes shown in figure 4.16 for this structure. Results are shown in figure 4.17. As one can see, the total absorption is relatively well computed even on the coarse mesh. In a second time, the same study has been conducted for the volumetric absorption in the GaAs. In figure 4.18, the convergence when increasing the polynomial order is shown. As one can see, the convergence of the volumetric absorption needs to use much higher polynomial order in order to converge, compared to the total absorption. Those results have a physical explanation. The volumetric absorption is a quantity computed inside the elements (via volumetric integration), which makes this observable a "near field" quantity, while on the other hand, the total absorption (as well as the transmission and the reflection) is computed on a surface (via surfacic integration) "far" from the nanotextures. This makes the variations in the near field much less impacting. In figure 4.19, we look at the convergence when refining the mesh. The results are well converged, even for the coarse mesh. The results on the very coarse mesh are even better than the results on the coarse mesh with a lower interpolation order. Remark : The volumetric absorption is, by nature, more expensive to evaluate than the total absorption. This is even accentuated with the need for higher polynomial order.

Results

All simulations have been run on 128 cores (16 Intel E5-2670 2.6 GHz processors, each with 8 cores). In figure 4.20, we plot A tot obtained for the DGTD method based on P 3 and P 4 interpolations using a coarse mesh. As one can see, A tot is quasi-converged for P 3 on every structures. The same observation is made for the total absorptions in each of the structures, only a few resonances located near the GaAs bandgap, i.e. after 800 nm, are still presenting discrepancies. Regarding physics, a considerable improvement is observed between Structure 1 and the 4 others (from 300 nm to 700 nm). This corresponds to the adding of the anti-reflective coating. Structure 2 presents some absorption gaps in the upper part of the spectrum, which are gradually filled by breaking the symmetry (Structure 3) and adding a z-invariance in the nanostructure (Structure 4). In figure 4.21, we plot A for each material, obtained for the DGTD method based on P 4 interpolation using a fine mesh. The important physical quantity here is the absorption in the GaAs. As one can see, the same observation as for the total absorption can be made. However, we seem to see a decrease in the absorption spectrum of the GaAs between Structure 3 and Structure 4, but an increase of the absorption of the Ag. This is not the desired behavior, and in this case, we see that the volumetric absorption allows us to determine it, which was hidden by the total absorption.

Finally in figure 4.22, we compare our result on the total absorption with the other selected numerical methods, i.e. RCWA (Reticolo), FDTD (Lumerical FDTD) and FEFD. This study is still ongoing, which explains the absence of some methods for a few structures. The results observed are generally close. Differences are mainly focused on the upper part of the spectrum, the one that raised convergence problems for our method. The biggest differences are found in structure 5, where the two time-domain methods (FDTD and DGTD) are very close in comparison to the FEM. For this structure, RCWA results were not available due to the complexity of the nanostructure.

Conical texturation

In this section, we study two solar cell designs, which have been proposed in recent years by researchers in the field. These solar cell devices exploit gratings with nanocones in order to increase their light trapping capability. The first one is a double grating presents in [WYL + 12]. The second one is a conical nanowire presents in [START_REF] Eyderman | Near perfect solar absorption in ultra-thin-film {GaAs} phtonic crystals[END_REF]. We use our high-order DGTD method, and the obtained results are also compared with results provided by the most widely used numerical methods in the field.

Double grating

There is a significant recent interest in designing ultra-thin crystalline silicon solar cells with an active layer thickness of a few micrometers. Efficient light absorption in such thin-films requires both broadband anti-reflection coatings and effective light trapping techniques, which often have different design considerations. In [WYL + 12], the authors show that by employing a double-sided grating design, one can separately optimize the geometries for anti-reflection thanks to the top conical grating and also light trapping purposes, allowing broadband light absorption enhancement, thanks to the bottom conical grating. In the present study, we try to reproduce this test case. This structure contains a crystalline silicon thin-film with two silicon nanocone gratings. The nanocones form two-dimensional square lattices on both the front and the back surfaces. The film is placed on a perfect electric conductor (PEC) mirror. The design is introduced in [WYL + 12] as an optimal configuration, for a fixed quantity of absorber (corresponding to an equivalent flat layer of 2 µm c-Si). In the latter, the optimization is done using a RCWA solver with the aim of maximizing the Jsc.

Mesh

One can see in figure 4.23 the design of the cell. The two gratings are placed on each side of a large silicon layer of 1736 nm. The heterogeneity of element sizes can be seen directly. This is mostly explained by the geometric singularity at the top of the cone. In order to reduce this phenomenon, we choose to smooth out the extremity of all cones. Since we change the topology of the device, we readjust the equivalent volume of Si by increasing the size of the central layer according to the amount of volume truncated. As long as the cones are meshed with linear elements, Gmsh can be used.

Numerical convergence

We first conduct a series of simulations for the optimized structure obtained in [WYL + 12]. For the top nanocones, the period is 500 nm, the base radius is 250 nm, and the height is 710 nm. For the bottom nanocones, the period is 1000 nm, the base radius is 475 nm, and the height is 330 nm. Two tetrahedral meshes have been constructed using the Gmsh software. Mesh M1 consists of 34,264 cells. The minimal, maximal and average length of a mesh edge in this mesh is respectively equal to 0.011 µm, 0.210 µm and 0.101 µm. Mesh M2 consists of 147,508 cells. The minimal, maximal, and average length of a mesh edge in this mesh is respectively equal to 0.0014 µm, 0.110 µm, and 0.064 µm. Results of simulations with mesh M1 and using the DGTD method with interpolation order 1 to 4 are shown in figure 4.24 (absorption spectra) and 4.27 (contour lines of the amplitude of the DFT of E at a particular wavelength). These results shed light on a potential difficulty that one may face when assessing the numerical convergence of simulations of light trapping in solar cells, i.e. the convergence of near field quantities.

Results

In many works, the quantity of interest (QoI) on which the numerical convergence analysis is based is the total absorption A tot . When the FDTD method is used, the convergence is assessed by refining the underlying cartesian grid discretizing the solar cell structure. In the case of the DGTD method, two parameters can be considered separately or together to monitor the accuracy of a numerical solution: the discretization step referred to as h and taken as the maximum of the cell-wise sizes; the order of the polynomial interpolation of the components of the electromagnetic field within each cell referred to as p (we assume here, a uniform interpolation order). In figure 4.24, we fix the mesh and increase the interpolation order. We observe that the DGTD-P 3 yields a converged solution using mesh M1 (DGTD-P 3 and DGTD-P 4 solutions are indistinguishable). However, the contour lines of the amplitude of the DFT of E in figure 4.27 clearly lead to a different picture as the DGTD-P 4 solution shows patterns that are not well resolved when using the DGTD-P 3 method. Now, we perform a single simulation with the DGTD-P 2 method using mesh M2, see figures 4.25 and 4.28. Whereas the obtained absorption spectrum is almost identical to the one obtained with the DGTD-P 3 method using mesh M1 (see figure 4.26), here again, we observe that the volumic field is more sensitive to the refinement of the mesh than the total absorption. We conclude this series of results by assessing the influence of the physical simulation time t by comparing the absorption spectra obtained with the DGTD-P 2 method using mesh M2, for different values of T (note that T = 2.0•10 -13 seconds has been used in all the previously discussed simulations). In order to have a more generic stopping criterion, we numerically evaluate the residual energy present in the system. This energy E is calculated every n time-steps, n ∈ N * . Once the incident field has been fully injected, the maximum energy is recorded. In the course of the simulation, the energy in the system decreases to zero. We define the threshold for the absorption computation, which corresponds to a certain percentage of max n>0 E. For example, we usually used 0.01% as a suitable starting parameter. It can be noted that decreasing this threshold is equivalent to increasing T . The obtained absorption spectra for different values of the threshold are shown in figure 4.29. The simulation physical and computational times corresponding to the different energy levels are shown in table 4.3 We also distinguish a very different behavior between the small wavelengths and the large wavelengths. Even for the largest threshold (1%), the absorption for the lowest wavelength part is already converged. For the upper part of the spectrum, we observe stronger resonances that are still appearing when the threshold is strongly reduced (0.001% or a 200 longer physical time). In figure 4.30, we show the difference between two consecutive absorption spectra, which further confirms the difference in the behavior of the two parts of the full spectrum. This can be explained by the large extent of the spectrum of interest (200 nm to 1200 nm) that complicates the setup of the simulation. Indeed, it is well known that a good solution to wave problems requires a certain number of points per wavelength. However, for a time-domain solver, multiple wavelengths are handled synchronously. Here, the ratio λ max /λ min is 6, which means that we should have roughly six times more degree of freedom for the lower part than for the upper one.

Therefore, we have taken a "frequency-domain" approach in order to reduce the overall computational time. We choose to split the simulation into two parts, one dealing with the small wavelengths, i.e. with a refined mesh, and high interpolation order, while the other can run on a coarser mesh with a lower interpolation order. Indeed, the physical simulation time required to capture the resonances of the upper part of the spectrum will remain the same, i.e. a threshold 1.50e-03 6.02e-14 4 h 04 mn 0.5% 3.00e-03 1.00e-13 6 h 45 mn P 2 -P 3 5% / 3.82e-13 7 h 55 mn 2% 5.22e-01 (4.28e-01) 5.23e-13 10 h 50 mn 1% 3.36e-01 (5.82e-02) 6.52e-13 13 h 30 mn 0.5% 1.79e-01 (3.71e-02) 7.93e-13 16 h 25 mn 0.1% 3.94e-01 (4.92e-02)

1.20e-12 24 h 55 mn 0.05% 1.76e-01 (1.01e-02) 1.38e-12 28 h 32 mn 0.01% 1.90e-01 (8.25e-03) 1.93e-12 40 h 03 mn 0.005% 1.11e-01 (1.90e-03) 2.31e-12 47 h 56 mn 0.001% 8.69e-02 (1.46e-03) 3.10e-12 64 h 12 mn Table 4.3 | Correspondance between the energy levels and the physical time for both spectrum area. The error in the parenthesis is computed on absorption after the moving average. lower than 0.005%. However, being able to use a much coarser mesh size, with a lower interpolation order will drastically decrease the computational load as well as increase the stable time-step. Also, in figures 4.31 and 4.32, we check the coupling between the absorption spectrum of the two simulations, as well as the concordance with the full run on M2. In figure 4.31, we superpose our obtained absorption spectrum with the single-pass absorption and the Yablonovitch limit, respectively given by A Single-pass = 1e -αd and A Yablonovitch = 1 -1 1+4n 2 αd . Finally, we have shown that the use of a time-domain method for "too much" broadband problems can pose complicated numerical problems. The use of a frequency-based strategy allows the simulation to be split into two parts, each of them being much cheaper in terms of computation cost. The overlapping between the absorption spectrum is almost perfect.

A tot (λ) E = 1% E = 0.5% E = 0.1% E = 0.05% E = 0.001% E = 0.001% (

Conical nanowire

In [START_REF] Eyderman | Near perfect solar absorption in ultra-thin-film {GaAs} phtonic crystals[END_REF] the authors study light trapping in Gallium Arsenide (GaAs) solar cells based on thin-film photonic crystals using a FDTD method. Indeed, light trapping in thin-film photonic crystals is capable of surpassing the performance of conventional thick solar cells. This is made possible by a combination of improved solar absorption and more effective carrier collection on length scales smaller than the carrier diffusion length. Thus, thin-film photonic crystal active layers provide a unique opportunity for both photonic and electronic management in solar cells. The authors present designs that enable a significant increase in solar absorption in ultra-thin layers of GaAs. In the wavelength range [400 nm, 860 nm], 90-99.5% solar absorption is demonstrated depending on the photonic crystal architecture used and the nature of the packaging. It is shown that using only 200 nm equivalent bulk thickness of GaAs, forming a conical-pore photonic crystal (lattice constant 550 nm, pore diameter 600 nm and pore depth 290 nm) packaged with SiO2 and deposited on a silver back-reflector, one can achieve 90% absorption of all available solar light in the wavelength range [400 nm, 860 nm]. We consider here the optimized setting studied by the authors, which is characterized by a 200 nm equivalent bulk thickness of GaAs (which represents a height of 4770 nm), a cone base radius of 100 nm, and a lattice constant of 500 nm. Such geometry is complicated to represent with structured mesh. In this part, we want to show the advantages of our DGTD method for simulating light absorption in such PV cell. One mesh of the structure is shown in figure 4.33.

|A i tot (λ) -A i+1 tot (λ)| E2 -E1 E4 -E3 E7 -E6 (

Numerical convergence

The device geometry presents real challenges. Indeed, the combination of a high height (h = 4.77 µm) and a small basis radius (r = 100 nm) for the GaAs cone requires the use of very small element at the top (h min = 2nm). Mainly due to meshing issues, we choose to simplify the extremity of the cone by truncating it. (see figure 4.34).

The convergence was tested on both fine and coarse meshes. The characteristics regarding the meshes are in table 4.4. We plot in figure 4.35 and 4.36 the total absorption spectrum obtained for polynomial orders going from 1 to 4 for the coarse and the fine mesh, respectively. As can be seen, for A tot , P 3 and P 4 simulations are really close.

Mesh

Tetrahedron h min (nm) h max (nm)

Results

We first discuss the numerical treatment of the optical characteristics of GaAs. This material exhibits a zero absorption above its bandgap which is around 870 nm. While this zero absorption feature is easily accounted for in the frequency-domain modeling setting by simply taking ε i (λ) = 0 for λ > 870 nm, it represents a more challenging task in the time-domain case. The generalized dispersion model presented in details in [START_REF] Viquerat | Fitting experimental dispersion data with a simulated annealing method for nano-optics applications[END_REF] involves several parameters. The process of fitting experimental optical data yields numerical values for these parameters through the resolution of an optimization problem. Simulated annealing is an appealing algorithm in this context, in particular, because it is a global optimization method. However, this algorithm works with several parameters, such as the number of neighbors used, the neighborhood size, or the temperatures, whose optimal values for a given problem and a target accuracy requirement in the solution of this problem, are in general difficult to infer. We have considered two parameter configurations (referred to as original and new in the following) and run the simulated annealing algorithm for fitting the real and imaginary parts of the electric permittivity of GaAs. The solution with our initial guess (option original), which works well for many other materials, is shown in figure 4.37. At first glance, the obtained result seems satisfactory, although, by zooming in the bandgap region, we clearly observe a mismatch in the fitting of the imaginary part. The left plot of figure 4.40 shows the result of a simulation with the DGTD-P 2 method in terms of the total absorption spectrum. We note that the time-domain simulation fails to accurately model the zero absorption above 870 nm, which fits with the mismatch of ε i . Indeed, the accuracy of the simulation of light absorption in the studied GaAs solar cell appears to be highly sensitive to the quality of the fitting, especially near the bandgap. We tried to tune the simulated annealing algorithm in order to improve the quality of the fitting, by adding some constraints on the bandgap region in the optimization part. One can see the clear improvement in figure 4.38. We can see a comparison in log-scale in figure 4.39, which clearly highlights the fitting issue. The simulation performed with the second fit confirms the importance of the modelization of the imaginary part of the permittivity, which produced almost three times more unphysical absorption (see figure 4.40). In this section we have shown the importance of the material model used, which has a huge influence on the final results. The bandgap's fit problem is a known problem in the time-domain field: in [START_REF] Eyderman | Near perfect solar absorption in ultra-thin-film {GaAs} phtonic crystals[END_REF], the absorption is plotted from 400 nm to 860 nm.

Thin-film solar cells

In this part, we present a joint work realized in interaction with physicists from IEK-5 Photovoltaic, Forschungszentrum Julich (FZJ). This work is part of the EoCoE project [START_REF] Lanteri | Rigorous modeling of light absorption in nanostructured materials using a parallel high order finite element time-domain technique[END_REF]. We study light trapping in a silicon-based thin-film solar cell setup that consists of several randomly textured layers (see figure 4.42). The focus is on amorphous and microcrystalline silicon (a-Si:H and µc-Si:H), which belong to the family of disordered semiconductors. The main characteristics of those materials is the structural disorder, which affects the optical and electronic properties significantly.

Geometrical models

Researchers from FZJ provide us with topographic data of each layer obtained from AFM. In order to generate a geometrical model, we apply the different strategies presented in section 4.3. Partial views of generated tetrahedral meshes are shown in figure 4.43. These meshes correspond to a subset of the full model, which has been used for preliminary tests and comparison of the different strategies for periodic boundaries. It is obtained by taking only 1 × 1 µm 2 subset of each layer data instead of the full 10 × 10 µm 2 points. The mesh in figure 4.44 is obtained for the full model using the Tukey window approach. Indeed, it is way too costly to use the alternative method for imposing periodicity as it yields a tetrahedral mesh with 5 million elements.

Material models

The optical properties of the different materials that constitute the considered solar cell structure have been fitted to the parameters of our generalized dispersion model, which was originally intended for metals. The obtained permittivity functions are plotted in figure 4.45. As can be seen, all materials are relatively well approximated in the range λ = [300, 1000] nm. Regarding glass, a constant permittivity ε r = 2.25 is used.

Numerical convergence

Numerical convergence has been assessed by considering the small submodel of the solar cell structure on the one hand, and the full-size model, on the other hand. The first mesh (M1) is used as a reference: the nanotextures are removed and the interfaces are flattened. Two tetrahedral meshes (M2 and M3) have been constructed on the realistic topography data, using the procedure described in section 4.3. The characteristics of the meshes used in this study are summarized in table 4.5. In this table, h min and h max respectively stand for the minimum and maximum length of a mesh edge. We plot in figure 4.46 to 4.48 the EQE, for both materials, obtained for polynomial orders ranging from 1 to 4 in the DGTD method. As can be seen, the DGTD-P 3 and DGTD-P 4 methods yield almost identical results when considering the submodel for the a-Si:H layer. Moreover, we observe a smoothing of the plots from M1 to M3. Indeed in figure 4.46, the EQE on M1 presents some resonances for both the µc-S:H (between 600 nm and 800 nm), and for a-Si:H (at 430 nm and 500 nm). We seem to be observing some residuals of this with the submodel in figure 4

Comparison with FDTD simulations

A comparison with results from FDTD simulations performed at IEK-5 is shown in figure 4.49. As can be seen, the results are in relatively good agreement, especially for the a-Si:H. The observed discrepancies can be attributed to two causes: (i) the use of a Cartesian grid in the FDTD method and (ii) the number of frequency points for the evaluation of the observable. Concerning the first point, the results of the FDTD simulations presented here have been obtained for a Cartesian grid that is refined in order to minimize the staircasing effect inherent to this type of discretization, while keeping the simulation time to a reasonable level. A more accurate comparison would require to refine the underlying Cartesian grid further, therefore, increasing the FDTD simulation time. Regarding the second point, the FDTD results have been obtained by considering a monochromatic plane wave for the incident field and sweeping over the target frequency range. Because of the computational overhead, the number of FDTD simulations has been kept low, as can be seen in figure 4.49. Shannon theorem for frequency-dependent observables Since we are considering a timedomain modeling setting, the computation of the observable quantities discussed previously requires to perform a discrete Fourier transform of the electromagnetic field on the fly. Therefore, in order to reduce the computational overhead of computing the volumetric absorption, we have exploited Shannon's theorem for sampling of frequency-dependent quantities. In fact, setting:

∆t obs = 1 2f max ,

Conclusion

In this chapter, we started by reviewing the three major numerical methods used in the PV field to deal with the optical part of PV problems. Then, we have presented the main quantities of interest used when accessing the quality of solar devices, precisely the external quantum efficiency (EQE) and the short circuit current (Jsc). The procedure to build a realistic mesh starting from a topography measured, generally obtained by atomic force microscopy (AFM), has been presented. Some geometry processing has been explained, in order to get usable meshes, mostly in terms of periodicity. Then we showed different realistic solar cells simulation, from several collaborations with physicists. We provided convergence analysis for all the cases, and we compared our results with the state of art solvers. Our methods obtained relevant results that are close to the other methods. In future work, it would be interesting to validate our results with experimental results. Also, a comparison in terms of computation time is mandatory in order to highlight which methods are the most profitable depending on the considered case. Optimization is an essential tool for making decisions. It focuses on finding the best solution among the set of all feasible solutions. It is used in many science and engineering fields. When coupled with numerical solvers, it allows prototyping at virtually no cost. For this reason, optimization is widely used in photovoltaics. Most of the time, the quantity to maximize is the absorption of solar light by the cell. One of the main challenges and research axis consists in increasing the light trapping in the absorber to exceed the best possible efficiency [START_REF] Garnett | Light Trapping in Silicon Nanowire Solar Cells[END_REF]. In [START_REF] Ganapati | Light Trapping Textures Designed by Electromagnetic Optimization for Sub-Wavelength Thick Solar Cells[END_REF], Ganapati et al use optimization in order to design nanoscale textures for light trapping in subwavelength thin-films. The optimization is based on an adjoint gradient method to search for a local optimum, and is coupled with FDTD optical simulations. The authors represent the surface texture as a truncated Fourier series, and the optimization parameters are the Fourier coefficients. In [START_REF] Dühring | Optimization of extraordinary optical absorption in plasmonic and dielectric structures[END_REF], Düring et al use topological optimization on the geometry of the surface structuring for both plasmonic and dielectric devices. The merit function is the extraordinary optical absorption (EOA) at one specific wavelength. They show that the topological optimization can propose non-intuitive surface designs. In [LYE + 17], Lee et al study the quasi-random features that can be imposed on an amorphous silicon surface by wrinkle lithography. They use Fourier-based inverse design and show an enhancement of light trapping capability by a factor of five over the [400 nm, 1200 nm] region. In [ACM], Anderson et al are optimizing the efficiency of thin-film PV solar cells. The optimization algorithm is based on a gradient-free optimization algorithm, which is well suited for maximizing an objective function over a high-dimensional parameter space, coupling with an HDG optoelectric solver.

Design optimization

Most of the optimization problems in science and engineering are black box optimization problems, and are characterized by an objective function f (x) ∈ R that does not have an analytic expression. The function f is depending on the parameter x, which is searched in a space called the admissible design space, denoted by W ⊂ R d . Evaluating f can be cheap, which allows a sampling of many points in W, e.g. , via random search, or it can be expensive, such as when it requires the simulation of a complex physical problem in three space dimensions, or when the data required to evaluate the objective function come from laboratory experiments. In such situations, one has to select an optimization technique that minimizes the number of evaluations of the objective function. Moreover, the problem at hand may exhibit several local minima, thus motivating the use of a technique that can deal efficiently with this issue and provide a global optimal solution. In this context, Bayesian Global Optimization (BGO) techniques are particularly relevant [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF].

In this chapter, we consider an optimization strategy based on an adaptive statistical learning approach. The method used is called Efficient Global Optimization (EGO). It was initially introduced in [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF], and it is now one of the most popular BGO techniques. This method relies on the construction of metamodels or surrogate models. Historically, linear regressions were used as metamodels, while today Gaussian Process (GP) models are used. It has been initially developed for optimization studies concerning fluid flow problems [START_REF] Duvigneau | Kriging-based optimization applied to flow control[END_REF][START_REF] Duvigneau | Comparison of turbulence closures for optimized active control[END_REF] and fluid-structure interaction problems [SHD + 17]. We briefly review the design and working principles of this optimization approach. GP based optimization methods have become very popular in recent years for the global optimization of complex systems characterized by high computational costs in the evaluation of the objective function. GP models have also been found especially interesting for optimization, particularly in the framework of the EGO method. One of the main ideas of the EGO is to rely on an acquisition function, usually, the Expected Improvement (EI) [START_REF] Mockus | Global Optimization and the Bayesian Approach[END_REF], to find new promising designs. During the iterative process, the new point obtained is tested in the objective function, and its result is used to update and improve the accuracy of the metamodel.

The chapter starts with a presentation of the EGO algorithm in section 5.1. The optimization notations are explained, and a simple example is provided to illustrate the method. Section 5.2 deals with the study of a textured PV glass. The main objective is to reduce the reddish reflection caused by the silicon while enhancing the total absorption. To do so, we optimize the pyramidal elementary pattern, which is replicated into all layers. Section 5.3 deals with a hot carrier solar cell. It is using a thin layer of indium gallium arsenide (InGaAs). The main idea is to work on the nanostructured back mirror to increase light trapping. This case is optimized using a parallel EGO (also named q-EGO), which allows to accelerate the iterative part of the EGO algorithm by evaluating multiple objective functions in parallel.

Efficient Global Optimization (EGO)

Bayesian Global Optimization (BGO) provides efficient and effective techniques based on Bayes Theorem to direct the search of a global optimization problem. To do so, it starts by building a probabilistic model of the objective function, called the metamodel or surrogate model. The specific BGO technique used in this work is called Efficient Global Optimization (EGO). EGO is used when the evaluation of the objective function is an expensive task, which is precisely the case in our context in which we resort to a full-wave electromagnetic solver in three space dimensions for each cost function evaluation. EGO consists of two steps. The first one is called the learning step. It relies on a set of values of the objective function at several points in the parameter space. Those values are used to build a surrogate model. It is usually called the Design of Experiments (DoE).

Design of Experiments

The initial metamodel is constructed from a design of experiment (DoE) denoted by X DoE = {x 1 , ..., x n } ∈ (W) n , where n ∈ N, n > 0 is the number of points in the DoE. Usually, W = R d . The set X DoE must be chosen depending on the complexity of both f and W as the accuracy of the resulting metamodel depends on it. Different possibilities are available to construct the design, such as uniform sampling, random sampling, orthogonal sampling, or the reference method, which is used in this work, the Latin Hypercube Sampling (LHS) [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. LHS works as follows: it positions each design in the DoE database such that it does not have common coordinates with the other previously positioned designs, as illustrated in figure 5.1. The main advantage of this strategy is the minimization of the interaction of the variables, as shown in [START_REF] Stein | Large sample properties of simulations using latin hypercube sampling[END_REF]. Then, one evaluates the objective function for each x i ∈ X DoE to get the set of data (X DoE , f (X DoE )). This allows to construct a surrogate model that aims at mimicking precisely the behavior of f while being computationally cheaper to evaluate.

Surrogate model

A popular surrogate model for Bayesian optimization is a Gaussian process (GP). GP modeling is a statistical method to approximate functions from a finite set of possibly noisy observations at ar-bitrary points. An initial surrogate is constructed using the dataset from the DoE. This metamodel is then used to predict at any promising design x t where t, t > n represents the index of the design. Then one evaluates f (x t ), which is added to the dataset X = {(X DoE , f (X DoE ), (x t , f (x t )}, and used again to update the surrogate model. GPs are used for both their cheap evaluations and their interesting statistical information. Indeed, one can predict the objective function at points x t thanks to the GP model mean, as well as an estimation of the error in terms of the GP model variance.

Acquisition functions

The new designs used to update the surrogate are obtained via an acquisition function, which must be defined in such a way that it allows both exploitation and exploration. Exploitation uses the predictions of the surrogate model to propose a new design maximizing or minimizing the objective function, while exploration proposes a new design in the area of the design space where the uncertainty of the model is high. More formally, the new design x t is defined by x t = argmax x∈W g(x|X 1:t-1 ) where g is the acquisition function and X 1:t-1 = {(x 1 , f (x 1 )), ..., (x t-1 , f (x t-1 ))} are the t -1 samples, computed from f so far with t > n. The most commonly used acquisition function, and also the one considered in this work, is the Expected Improvement (EI). To give an idea of what is the expected improvement, we place ourselves in the case where we have evaluated t -1 designs. The goal is to maximize f on W, and our current largest evaluation is

f + 1:t-1 = max x∈X 1:t-1 f (x).
We suppose that we have one additional evaluation x t to perform. After this evaluation, the new optimum will either be f +

1:t-1 if f + 1:t-1 > f (x t ), or f (x t ) if f + 1:t-1 < f (x t ).
One can quantify the improvement over this additional evaluation as max(f (x)f + 1:t-1 , 0). As f (x) is unknown until its evaluation, one can take its prediction from the surrogate model and so choose x t . More formally, the EI is defined as

EI(x) = E max(f (x) -f + 1:t-1 ), 0), (5.1)
which can be evaluated analytically under the GP model as

EI(x) = (µ(x) -f + 1:t-1 )Φ(Z) + σ(x)φ(Z) if σ(x) > 0 0, if σ(x) = 0, (5.2)
where

Z =      µ(x) -f + 1:t-1 -ξ σ(x) if σ(x) > 0, 0 if σ(x) = 0,
(5.3) with µ(x) and σ(x) are respectively the mean and the standard deviation of the GP prediction at x. The functions Φ and φ are respectively, the cumulative distribution function (CDF) and the probability density function (PDF) of the standard normal distribution. The first term in equation (5.2) is the exploitation term, and the second term is the exploration term.

EGO workflow and illustrative example

The DoE phase is carried out to construct an initial model. Then, an iterative and adaptive strategy is used to enrich the model and locate the most interesting values of the cost functional.

In this perspective, an internal optimization problem is solved to determine a set of new promising designs to be evaluated. The new cost functional evaluations are then employed to update the model. This iterative enrichment step is conducted until a stopping criterion is verified. It can be an a priori number of iterations N , a budget time constraint e.g. the optimization should take less than 10 hours, or even a convergence threshold, e.g. based on the merit function as EI < . We illustrate in figure 5.2 the main steps of EGO. The operation of the EGO algorithm is now illustrated using a 1D example, which consists of maximizing an analytic objective function given by

f (x) = -cos 3x cos 2(x + 1) - x 2 3 , x ∈ [-2.2, 2.2].
(5.4)

Let us choose a DoE composed of two points, i.e. X DoE = {-1.6, 1.6}. We represent four steps of the EGO procedure in figure 5.3. On the left column, one can see the objective function (dashed blue line), the measures (black dots), and the surrogate model represented by its mean (orange curve) and its confidence interval (light orange surface). In the right column, one can see the expected improvement in orange along with x + . We represent in figure 5. 4 (a) the distance between two consecutive points predicted from EI, and in figure 5.4 (b), the current maximum value. It is interesting to see how the surrogate is evolving: at first, it seeks and searches in uncertain regions (from iteration 1 to 5), and then it starts focusing the refinement close to the actual maximum. This is one of the very interesting property of EGO, which allows to reach global optimums without being stuck in a local one. In the left we can see the distance between two consecutive x + . The iterations 1 to 5 correspond to the exploration part. In the right we plot the current maximum.

In the following, we apply the EGO method to the design of nanostructured solar cells for optimizing their light trapping properties. We consider two solar cell configurations corresponding to two different forms of nanostructure. Moreover, we consider EGO implementations from two different software frameworks: the FAMOSA software suite and the DiceOptim package for R.

FAMOSA The Full and Adaptive Multi-Level Optimum-Shape Algorithm 1 (FAMOSA) is an optimization platform developed at Inria Sophia Antipolis -Méditerranée by the Acumes project team. It is devoted to multidisciplinary design optimization. This software must be coupled to a solver that provides a merit function. FAMOSA can be coupled to any solver and, therefore, can be used for several purposes, ranging from numerical to multidisciplinary design optimization. It has been developed initially for aerodynamic design.

DiceOptim DiceKriging and DiceOptim2 are R packages available on the official R(CRAN) repository. They have been developed in the frame of the DICE (Deep Inside Computer Experiments) consortium. DICE gathered major French companies and public institutions, having a high research and development interest in computer experiments. The main objective was to put in common industrial and academic problems in order to foster research and transfer in the field of design and analysis of computer experiments.

Surface texturing for enhanced color rendering

This study has been done in collaboration with Sunpartner Technologies, a company designing and manufacturing PV cells embedded in mainstream technologies. In this context, we consider a silicon based PV glass. This glass is actually a multilayered structure in which each layer of material contributes to the efficiency of the whole device. This effectiveness is assessed as the capacity of the structure to reflect specific colors of the visible spectrum in variable ambient lighting conditions, and to absorb further colors that affect user comfort. A difficulty that Sunpartner Technologies is currently facing in the development of its PV glass is the inadequate absorption of red tints, as shown in figure 5.5 left, which cause visible red reflections. Those reflections come from the absorption coefficient of the silicon. Indeed, it is significantly higher in the blue part than in the red part, which results in a significant penetration depth of nearly fifty microns. Minimizing the reflection of the red part of the spectrum and, at the same time, improving the absorption properties of the device is achievable from two main strategies. The first one is to increase the thickness of the silicon layer to reach its penetration depth. The second one consists of embedding the cell with anti-reflection coating and light trapping textures, which are reducing reflections, and increasing light trapping, respectively. From an industrial point of view, the first option may not necessarily be economically viable. The alternative option is to use thin layers. However, when the PV cell structure consists of an association of thin-films of some hundreds of nanometer for the thickest, there may be creation of colors by constructive interference during the light propagation in the multilayer structure.

The surface texturing is the most promising solution to deal with both issues: it can limit the effects of these interference phenomena on one hand, and increase the poor absorption of predominantly red colors on the other hand. Figure 5.5 right illustrates three strategies for texturing a typical structure, representative of a PV glass. This structure is made of three materials (glass, silicon (Si) and Transparent conductive oxid (TCO)), possibly with an anti-reflective gel coating. In the following, we propose to study the impact of textures by modeling and numerically simulating the propagation of solar light in this type of PV glass.

Objectives of the study

In this study, we limit ourselves to an idealized PV glass structure based on an average topography of the nanotexturing of the different layers, as shown in figure 5.6 for a classical solar cell stack. For the PV glass structure considered in this study, Sunpartner Technologies has provided an average topography representative of actual built structures, with three examples shown in figure 5.7. The use of a PV glass structure based on an average topography allows us to apply a hypothesis of periodicity more conveniently, as shown in chapter 4, reducing the size of geometric models and, thus, the computation cost of the full-wave solver. In section 5.2.4, a parametric study is first carried out in order to obtain a preliminary picture of the role of the texturing on the absorption of red tints and the creation of colors by constructive interference in the layers. Then, in section 5.2.5, an inverse design analysis is conducted using EGO in order to obtain the optimal nanotexturing for a given objective.

Cell composition

We consider a PV structure made of 7 layers (see figure 5.8). Layers characteristics are summarized in table 5.1. We can see that this cell is alterning very thick layers of hundreds of nanometers with very thin layers of few nanometers.

Construction of geometrical models

In this study, instead of considering realistic texturing of the PV layers as depicted in figure 5.7, we adopt a homogenized texturing. An analysis of the texturing of several tens of cells allowed to choose the pattern used, which is a periodic pyramidal PV pattern. Pyramid parameters such as heights and bases, are deduced from surface roughness data given by Atomic Force Microscopy (AFM) (defined in chapter 4). A geometrical model of this structure is shown in figure 5.9. Such simplified textures are useful for numerical simulations because they naturally allow the use of PBC, thus reducing the size of the computational domain drastically. To do so, we need to construct the elementary pattern of the structure, as shown in figure 5 even more the size of the computational domain to one-quarter of the elementary pattern if normal incidence is assumed, see figure 5.10(b).

Remark : Those geometrical simplifications are possible since we are not considering the realistic topography, as in figure 5.7. However, it would be interesting to validate the reliability of the simplified model with the realistic one, which will be the subject of a future study.

Meshes are obtained using the Gmsh3 CAD modeler and mesh generator. The elementary pattern that we consider in this study is modeled as the union of five pyramids.

Remark : Such a geometrical model is not trivial to create, especially when the Computer-Aided Design (CAD) has to be modified and rebuilt in the framework of an optimization process. Indeed, in this context, the parameters are supposed to change almost freely. However here, if pyramid height parameters are bigger than the upper layer starting altitude, pyramids can intersect pyramids from the upper one, changing most of the underlying geometrical topology. Such issues are generally hard to handle and often require manual modifications.

Influence of the geometry of the pyramids

Based on the results of a preliminary study made by Sunpartner Technologies, which was concerned with the measurement of texture roughness of the PV glass samples, a particular texture has been selected for the purpose of this study (see figure 5.11). The characteristics of this texture are:

• Ra (nm): 47.8 ± 3.7 (average roughness);

• Rq (nm): 63.6 ± 8.4 (average quadratic roughness);

• Sal (nm): 165.1 (auto-correlation length).

The parameter Sal corresponds to the distance over the surface such that two points have minimal correlation. This value is used to set the periodicity of the simplified geometry. The parameters Ra and Rq correspond to the roughness of the texture. Those three parameters are used to determine the general texture specification (see table 5.2). In term of numerical and geometrical model, we use a layer of Perfect Match Layer (PML) at the bottom of the cell to remove the parasitic Material h (nm) r (nm) AZO 47.8 165.1 reflection. The incident light is propagating from the bottom to the top. In case of the periodic set-up, PBCs are used in both xy directions. For the quarter domain, PEC and PMC are used, with PEC in the polarization direction. The top of the cell corresponds to the Al back reflector, which is completed with a PEC boundary.

We first demonstrate the influence of the interpolation order in the DGTD method on the convergence of the absorption spectrum on a chosen geometry (h and r are fixed). For this purpose, we use a coarse tetrahedral mesh with only 1001 cells, and we consider the geometrical characteristics of pyramids given in table 5.2. The absorption spectra resulting from simulations with the DGTD-P 2 and DGTD-P 3 methods are shown in figure 5.12. Simulations have been performed on a multi-core server equipped with Intel Xeon Gold 6154 CPU@3.0 GHz. Using four cores the simulation time is respectively equal to 385 sec (DGTD-P 2 method) and 1181 sec (DGTD-P 3 method). We see that the solution resulting from the DGTD-P 2 run is already sufficiently converged everywhere. An absorption above 60% is observed from 300 nm to 650 nm. The problem concerning the red part is clearly illustrated here: the absorption drops below 40% around 700 nm to 750 nm. Nevertheless, we choose to adopt the DGTD-P 3 method in the sequel in order to be more confident with the convergence of our results.

In a second step, we perform a preliminary parametric study in order to assess the influence of the geometrical parameters h and r. Various configurations have been tested and are summarized in table 5.3. The corresponding absorption spectra are plotted in figure 5.13, 5.14 and 5.15. From figure 5.13, we restrict the modification to the pyramid radius. We observe a small improvement for the absorption around 650 nm. It even fills a gap at 670 nm. However, it did not bring any improvement in the upper part of the spectrum. In figure 5.14, we do the counterpart by changing only the height of the pyramid. We see some interesting improvement on the upper region of the spectrum with a slightly decay for the lower part. For the results in figure 5.15, we take two other couples of (h, r), both higher than the starting ones. Combining the two options seems to further increase the absorption in the upper red part, with again, a small decay for the resonances before 600 nm. From those statements, we are confident with the obtained results which confirm the effectiveness of the optimization strategy.

h (nm) r (nm) # elem h max /h
A tot (λ)
h=47.8 nm, r=165.1 nm h=47.8 nm, r=82.5 nm h=47.8 nm, r=214.6 nm

Optimization of the geometry of the pyramids

Now that we have a preliminary picture of the influence of the two parameters h (nm) and r (nm) on the absorption profile, we realize an inverse design analysis. In this part, we make use of the EGO method from the Famosa library. Starting from the geometrical characteristics of the average pyramidal grating given in table 5.2, we investigate here the possibility of improving further the performance of the PV glass structure by using a numerical optimization strategy. In order to reduce the red reflections while increasing the total absorption, we choose to minimize the total reflectance over the spectrum of interest, i.e. λ ∈ [300 nm, 1000 nm], i.e. More precisely, we minimize the Euclidean norm of 500 equally spaced discrete values of the absorption in this wavelength range. The inverse design analysis has been conducted on a server equipped with 8 Intel Xeon Gold 6154 CPUs running at 3.00 GHz. The DGTD-P 3 method has been used for all the simulations. Given the number of parameters (two), and the size of the spaces for these two parameters, the DoE database is made of 20 points. For the same reasons, the EGO algorithm is set up for 20 iterations. Three optimization scenarios are considered, each with its own accessible design spaces. We have the scenario V1 with h ∈ [40 nm , 80 nm] and r ∈ [120 nm , 200 nm], scenario V2 with h ∈ [40 nm , 120 nm] and r ∈ [120 nm , 300 nm] and scenario V3 with h ∈ [40 nm , 180 nm] and r ∈ [120 nm , 300 nm]. The parameter ranges are considered as constraints in the EGO algorithm. The results of these three optimization scenarios are shown in figure 5.16 (a) to 5.16 (c). For scenario V3, the simulation time for the DoE phase is 5 h 24 mn, while for the EGO phase, it is equal to 6 h 36 mn.

Remark : The DoE phase is computed sequentially, which is not the optimal way. Indeed, one can compute all the DoE points in parallel, and then start the iterative part of the EGO. However, in our case, since the full optimization (DoE + iterations) is submitted entirely on 8 cores, it is easier to have the same behavior for the two parts, i.e. a sequential optimization process, while fullwave DGTD simulations are distributed on 8 cores. For V2 and V3, we observe a significant improvement for the absorption of the red color (≈ 700 nm). In the range [700 nm, 1000 nm], an improvement of almost 60% is noticed. In addition, there are no more wavelengths for which the absorption is lower than 50%. The contour lines of the amplitude of the DFT of E at two particular wavelengths are shown in figure 5.17. We clearly see what range of frequency each layer is targeting.

In conclusion of this study, using the EGO method from the FAMOSA library has allowed us to find better configurations both to reduce the reflection in the [700 nm, 900 nm] range and increase the overall absorption of the glass. Indeed, our figure of merit was not actually focusing on the reduction of the reflection. It could be interesting to consider this problem as a multi-objective optimization, with the second objective being the reduction of R in the range [700 nm, 800 nm]. Besides, a validation of the obtained optimal configuration is necessary. As we did in chapter 4 with the thin-film tandem solar cell, we could consider the exact topography for each layer, which would drastically increase the size and the complexity of the resulting geometric model. However, this comparison is necessary in order to show the reliability of the pyramidal-based simplified model.

Hot carrier pyramidal nanostructured solar cell

In the field of PV solar cells, many strategies are investigated for overcoming the Shockley-Queisser limit (defined in chapter 4). Multijunction solar cells are intensively studied for their capability of each individual junction to convert a fraction of the solar spectrum and give the record power conversion efficiency of 46%. However, producing cells reaching such values is extremely expensive, mainly due to the use of III-V materials. To reduce this cost, many studies have been performed on tandem solar cells. Hot carrier solar cell is an advanced concept which aims at overcoming the Shockley-Queisser limit with a single junction device where dissipation of the carrier energy via heat is avoided.

In this section, we consider a hot carrier solar cell proposed by researchers at C2N. The PV nanostructure considered is composed of several layers. From top to bottom: we have a 90 nm layer of Magnesium Fluoride (MgF2), then a 40 nm Zinc Sulfide (ZnS) layer. Just below, we have a 10 nm layer of constant refractive index n = 3.3, then the absorber layer of Indium Gallium Arsenide (InGaAs) of thickness 50 nm, another 10 nm layer of constant index. Below, we found pyramidal nanostructure, with a constant index n = 1.5. These pyramids are in a matrix of constant refractive index n = 3.3, and finally, we have a semi-infinite Ag layer at the bottom. A sketch of the cell structure is shown in figure 5.18.

Particularly in this structure, the thickness of the absorber plays a decisive role in the mechanism of the solar cell. Indeed, in order to avoid the thermalization rate, which reduces the cell efficiency, and improve the fast carrier collection rate, which increases the cell efficiency, it is necessary to use an ultra-thin absorber layer (less than 100 nm thick). In comparison, the previous cells (section 5.2) was composed by, at least, 400 nm absorber. However, such a thin layer drastically decreases the absorption of the cell. The main direction for overcoming this issue is to rely on an efficient light trapping strategy. In this work, we consider a back mirror nanotextured made of pyramids. This geometric model was chosen for its asymmetry, which promotes the light trapping. Indeed, in [START_REF] Cai | Effect of symmetry in periodic nanostructures on light trapping in thin film solar cells[END_REF], it is proven that both the level of absorption and the number of resonances increase by symmetry breaking. The goal here is to find the perfect pyramidal configuration, such that it increases the Jsc of the solar cell. To do so, we optimize the design of pyramid's grating, which is characterized by three parameters: the height of the pyramids h, the diameter d, and the gap between two elementary patterns g. The parameters are represented in figure 5. [START_REF][END_REF]. In comparison with the previous study, here the asymmetry of the pattern requires the use of PBC.

Material models

Concerning the optical characteristics of the constituting materials, for simplicity, we only consider the absorption in the InGaAs and the parasitic absorption in Ag. The absorption in the MgF 2 and the ZnS are neglected. All the fits are represented in figure 5.20. Concerning the imaginary part of the permittivity function of the InGaAs, we can see in figure 5.21, which represents the imaginary part of the permittivity in log-scale, that the representation of the bangap is problematic. As the fitting of silicium in the previous chapter, we tried to work on the simulated annealing algorithm to improve the accuracy near the bandgap. However for this one, it was not possible to improve it. It is certainly due to the very large spectrum of interest [300 nm,1600 nm], which is harder to take completely into account.

Comparison with FDTD results

Thanks to a preliminary multi-parametric study using a FDTD fullwave solver (Lumerical FDTD4 ), our physicist partners have proposed two potential designs: (D1) h = 390 nm, d = 690 nm and g = 0 nm; (D2) h = 390 nm, d = 690 nm and g = 8 nm. D1 and D2 correspond to the best designs obtained via the multi-parametric study.

In a first step, we reproduce the FDTD results for configurations D1 and D2 with our DGTD method. We performed a convergence analysis for both configurations (see figure 5.22).

The overall absorption profile obtained with both methods matches. However, the FDTD solution is closer to the DGTD-P 2 simulation, which is not converged. The difference between the short circuit current density (Jsc) computed with DGTD-P 2 and DGTD-P 3 is 12% (from 39.8 mA / cm 2 to 34.3 mA / cm 2 ). The other computed Jsc are in the table 5.4. We also plot two theoretical limits, the single-pass absorption and the Yablonovitch limit (defined in chapter 2) in figure 5.23.

Parallel EGO

From now on, we only consider the PV cell geometry without a gap. One limitation of the EGO from the computational efficiency perspective is the sequential evaluation of f in the optimization step, caused by the single points given by the EI figure of merit. Ginsbourger et al proposed in [START_REF] Ginsbourger | A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes[END_REF] a generalized acquisition function well-suited for parallel computing. This generalization is called multi-points expected improvement, or q-EI, as it generates q design points at each iteration of the EGO algorithm. This acquisition function allows us to perform the most expensive steps of the EGO, which is the objective function evaluation in a parallel way. Assuming that a minimum number of designs are necessary to converge towards the global optimum, the use of q-EI allows us to treat q points in parallel, and thus, to update the Gaussian model with q observations rather than one. Doing so, we can speedup the enrichment phases by evaluating multiple good candidates at once, and thus significantly improve the model in one iteration. An EGO diagram with q-EI is illustrated in figure 5.24. We have done the optimization of the pyramidal grating again, this time using the q-EI, also available in the DiceOptim library. The difficulty in this situation is the wrapper around the DGTD solver runs. Indeed, launching several instances of the DGTD solver in parallel, knowing that the solver is itself parallel, requires a careful management and distribution of computational resources.

The following test was performed on a cluster of Intel Xeon Gold 6154 CPUs running at 3.00 GHz. We start from a DoE database with 16 points, and we choose to predict and evaluate four points per iteration, over five iterations. As a comparison point, we also perform a sequential optimization in order to compare the obtained results for both methods. In order to also compare the optimization time, we take the same computation power per DGTD solver run (32 cores).

The results of different steps of the parallel EGO are shown in figure 5.25. Those maps represent the prediction of the metamodel at different updates of the model. For this structure, the maximum Jsc obtained for the optimized parameters d = 943 nm and h = 877 nm is 42.43 mA / cm 2 compared to 34.32 mA / cm 2 for the initial configuration. It corresponds to an improvement of 23%, as seen in the absorption spectra in figure 5.26. We also plot the contour lines of the amplitude of the DFT of E at wavelength λ = 1060 nm in figure 5.27. However, we note a nice improvement in the higher part of the spectra, which corresponds to the lack of accuracy of our InGaAs dispersion model.

Find x + = argmax x∈W EI(x) Solve f (x + ) Stop criteria Update GP model Find x + = argmax x∈W EI(x) Solve f (x + ) Stop criteria Update GP model Find x + = argmax x∈W EI(x) Solve f (x + ) Stop criteria Update GP model Choice of DoE Evaluate f (X) Construct initial surrogate model Find x + = argmax x∈W EI(x) Evaluate f (x + ) Stopping criterion Update surrogate model Stop • Optimize f (x) for x ∈ W • f is an expensive black-box

Performance

A comparison has been made between parallel and sequential optimization workflows. As we use four times more resources, we expect to have a computation time four times less important, assuming that there is no extra cost on the metamodelization side. We obtain a CPU time of 21 hours for the parallel version and 62 hours for the sequential one, which corresponds to a speedup of almost 3. We can explain this lower acceleration by a load balancing issue between concurrent evaluations of the solver. Indeed, designs produced by q-EI at each iteration can show substantial differences, resulting in entirely different computational domain sizes. The space of accessible designs being large, there can be a factor of 1000 between the volume of two structures proposed, i.e. (h pyr , d pyr ) = (100, 100) and (h pyr , d pyr ) = (1000, 1000). We show in figure 5.28 two designs evaluated during the last iteration. Design (a) corresponds to an exploration where

Case

Solver Jsc (mA / cm there has not been much research, while design (b) exploits the metamodel by searching in a potentially interesting area. The two simulation times associated with these problems are 33 sec and 11783 sec, respectively. An intelligent way to manage these problems would be to allocate computational resources not uniformly across simulations, but rather according to the complexity of the corresponding problem. Another way would be to use an asynchronous EGO algorithm, which updates the metamodel each time an evaluation of the cost function is completed.

Conclusion

In this chapter, we have presented an optimization method particularly well adapted to the typical objective function obtained by solving long simulations of DIOGENeS-DGTD. We processed two realistic PV cases with two different software environments (Famosa and DiceOptim). The first case consisted of a PV glass. Such technology requires a good efficiency with good aesthetic properties. That is why it is crucial to reduce the parasitic and reddish reflections. Using the Euclidean standard as a merit function, we observed a significant improvement in the red part of the absorption spectrum. Furthermore, the general aspect of the absorption was smoother, which corresponds to an attenuation of parasitic reflections. The second case consisted of a hot carrier solar cell. To increase the absorption of the cell with respect to the Jsc, a nanotexturing tetrahedral pyramid grating was placed at the bottom of the cell for light trapping purposes. A parallel EGO method was then used to optimize the grating. This parallel method, contrary to the sequential version presented, allowed updating the metamodel through multiple evaluations at each iteration. A sub-optimal scalability was observed, which can be explained by the couples of design given by the acquisition function. One way to overcome this can be to use an asynchronous EGO algorithm or to modify the distribution of computing resources by performing a prior analysis of the complexity of the design to be simulated, for example, by looking at the number of cells for the given mesh. A recurrent concern with the numerical simulation of realistic 3D problems is the time to solution, necessary to obtain sufficiently accurate results. In many cases, the time required to obtain these solutions is not problematic. In other cases, a budget of time and computing resources can be assigned to the simulation, thereby limiting the simulation parameters (mesh size, interpolation order, physical simulation time) in order to match the budget. For example, optimization problems (as those discussed in chapter 5) need multiple full-wave solver evaluations, therefore, increasing further the wasted time in case of poor solver performances.

High performance computing aspects

The aspect related to the study of the performances of numerical software, as well as its modification to improve its computational efficiency, is part of what is called high performance computing (HPC). HPC stands at the crossroads of different scientific disciplines and skillsets. For example, the discretization of a PDE system by FEM yields the resolution of linear systems. Developing efficient and scalable linear system solvers requires ingredients from numerical analysis, linear algebra, and computer science.

In general, a numerical solver is evaluated by its efficiency. There are two types of efficiencies: the first one is whether the obtained solution is accurate enough, which is most of the time, the principal (and only) concern for scientists; the second one is the computation time required to obtain the solution, which can be utterly prohibitive in some cases. Using poorly adapted methods and very low-performance software can make some problems almost unsolvable.

Parallel computing is nowadays ubiquitous in large-scale numerical simulations of complex physical problems. As its name suggests, a parallel software is able to run on many processors in parallel. This makes it possible to use the largest and most powerful computers, often called supercomputers. These are usually based on a cluster of computing nodes, which are basically multiple computers connected together. Those machines can make their Central Processing Units (CPUs) work on the same problem at the same time. Multiple forms of parallelism exist in a modern supercomputer, which materialize at different levels of the software. The latter are implemented by different parallelism protocols, the most renown being MPI (for Message Passing Interface), which was intended to distributed memory systems, and OpenMP for multithreading and shared memory programming. One other benefit of MPI programming is to remove the impossibility to run an application on one system because of the memory limitation by splitting this memory across several CPUs. For instance, this can be very helpful when large linear systems of equations have to be solved.

Parallel softwares introduce a new criterion to measure software performances. Indeed, a good parallelization should, in theory, make it possible to reduce the computation time proportionally to the number of tasks used. The main issue that can appear is referred to as load imbalance. A proper load balancing consists of distributing approximately equal amounts of work among tasks so that all tasks are kept busy at the same moment. It is primordial for parallel programs when synchronization is needed. Of course, if all tasks need to communicate at one point, this operation will be performed when every task is ready, i.e. when the slowest one has finished its assigned work. As a result, having a perfect speedup is not straightforward, and some implementation adjustments are required to balance the computational load between each CPU. Besides, we should not neglect the overhead introduced with the parallelism, e.g. the communications introduced by MPI. Two types of parallelism regions can be defined depending on the ratio between communication operations and local computation: the fine grain parallelism, which corresponds to piece of codes doing small amounts of computation between communication operations, and the coarse grain parallelism, which represents part of code where the computational work is significant compared to communications.

In this chapter, we focus on the performance of the DGTD solver developed in this thesis, called DIOGENeS-DGTD, which is a component of the DIOGENeS software suite1 . In section 6.1, we study the parallel performances of the coarse grain parallelization of DIOGENeS-DGTD, which is based on SPMD strategy combining a partitioning of the underlying mesh and a message passing programming with MPI. We highlight some load balancing issues in the DGTD time loop in the case of a realistic solar cell simulation. We also show that those load balancing issues are even more prominent for simulations on hybrid meshes. We propose a strategy to reduce these load balancing issues, which consists in an effective partitioning method based on an iterative algorithm over simulations. This allows us to find the most expensive area without any a priori knowledge. In the second part of this chapter, we present a modification of the DIOGENeS-DGTD solver to include hybrid parallelization in order to relax the constraint on the mesh partitioning when using a large number of tasks.

Coarse grain parallelism

In the vast majority of realistic studies of light trapping in solar cells, numerical simulations cannot be performed in a reasonable time sequentially. This is due to the geometric complexity yielding computational meshes with many cells on the one hand, and to a large number of resonant modes in the structure requiring a very long physical time in order for their accurate resolution on the other hand. One efficient parallelism paradigm for discontinuous Galerkin methods is built on top of graph theory. Indeed, in order to distribute the problem on several computational nodes, the involved data structures, which are typically restrained to the mesh T h for DG methods, need to be dispatched across the available resources. This is possible thanks to the partitioning of the mesh. DG discretization is known for being particularly well suited for distributed memory parallelism. This is due to the compactness of the discretization, which produces a limited amount of information exchanged between mesh partitions. Of course, the only information to be communicated between neighboring partitions is linked to the d.o.f. of the boundary cells.

Mesh partitioning

We start with some notations. As in all the previous chapters, we denote the considered mesh by T h . Then, this mesh is split into n, ∈ N, n ≥ 1 submeshes, denoted by T i h , i ∈ {1, ..., n}. In the parallel computing context, n also represents the number of cores. The number of neighbors of each T i h , i ∈ {1, ..., n} corresponds to the amount of inter-partition communications required to share information from one submesh to another. So it makes sense to try to get partitions that are almost equally weighted (i.e. with the same number of cells) while minimizing the number of neighbors. It is also important to minimize the size of the neighborhood region (which represents the memory size of the communication). Of course, an MPI parallelization is, by nature, artificially increasing the amount of d.o.f. of the problem due to more ghost cells at the inter subdomain boundaries. So minimizing the number of partition cuts is essential in order to have good performances. The main metric of parallel efficiency is the speedup. It defines how the computational time is reduced with the number of CPUs, for a fixed total problem size, i.e. the ratio t on n cores t on 1 core . To do so, we start by representing T h as a graph, with its elements K being the graph nodes while the faces of K are graph edges. This graph needs to be optimally partitioned, that is with an equal number of nodes within each subgraph and minimum number of partition cuts. In order to do so, we interface DIOGENeS with the MeTiS library [START_REF] Abou-Rjeili | Multilevel algorithms for partitioning power-law graphs[END_REF], which is a collection of programs that can be used for partitioning unstructured graphs both on serial as well as on parallel computers. The partitioning workflow is represented in figure 6.1.

Parallel load balancing

One major issue faced while dealing with realistic problems is the complexity and heterogeneity of the domain. Additionally, one simulation may involve various modeling and numerical features (e.g. hybrid cells, locally adapted interpolation order, PML, dispersive materials) and may also need additional computations (quantity of interests, injection of sources). Each of these ingredients represents a specific computational cost, both in terms of memory consumption and in arithmetic operations. Considering the aforementioned possible peculiarities, the risk of having some load unbalancing in the computation is high. Different strategies are possible to handle the load imbalance. The first one considered here consists of adapting the partitioning strategy. In this section, we profile DIOGENeS-DGTD on several test cases presented in chapter 4. Then some partitioning algorithms are used to obtain partitions of this dual graph (c). From here, the submeshes are recovered.

We first illustrate a typical load balancing problem in our applications. In order to approach a realistic situation, we enable the computation of the reflection coefficient using the Shannon-Nyquist sampling theorem. The simulation is run for 100 iterations using a DGTD method with a P 4 interpolation. The mesh is made of 2976 elements. It is partitioned into eight subdomains using a uniform distribution of weights associated with the nodes of the graph. The partitioned mesh is represented in figure 6.2. The subdomains have from 1 to 4 neighbors, and the number of cells deviation from the bigger submesh to the smaller submesh is 1.02. Those statistics are convincing at first.

The profiling of DIOGENeS-DGTD is done using two tools developed at the Barcelona Supercomputing Center (BSC). The first one, Extrae2 , generates traces of the execution of the program. Tracing allows us to collect low-level information, directly on the execution of functions and also communications and idle times. The other tool is Paraver3 , which is used to analyze these traces.

We focus on the time loop, which is the outer iterative part of the solver. In a realistic simulation, this loop represents 90% of the running time. We represent on figure 6.3 the 14 steps computational time using hybrid meshes in sequential execution mode. Here we consider a smaller example from chapter 4. We consider a solar cell with the back and front texturing. In figure 6.4, we show the two meshes used to compare the scaling of DIOGENeS-DGTD on tetrahedral and on hybrid meshes. The two meshes are constituted of 33200 cells and 23000 cells, respectively.

The partitions are here built without any constraint. The parallel speedups plotted in figure 6.5 are satisfying at the beginning, with a downfall for the hybrid case. We even see that at some point, the DGTD solver on hybrid meshes is less effective than on tetrahedral meshes. Two main reasons are explaining this behavior. First, the tetrahedral mesh involves twice as many degrees of freedom as the hybrid mesh, thus improving the ratio of computing time to communication time. Secondly, communications and computations caused by hybrid faces are considerably more expensive. If the partitioning is done in this region, those heavy computations are distributed over the cores, causing a significant overhead.

It is relevant that the hybrid interfaces introduce extra computations. As seen in chapter 3, the elementary matrices must be computed on all the hybrid cells (a hybrid cell has at least one hybrid face; a hybrid face is shared by two cells from a different type). In the case of some hybrid artificial interface, extra information needs to be communicated, resulting in an extended buffer. In figure 6.6, we present performance figures for simulations based on constrained partitions, which have been obtained by setting a high penalty for the hybrid faces weights. In this way, configurations with cutting on non-conforming faces are prevented. We observe a clear improvement of the scaling, which keeps more advantageous the use of a hybrid mesh than the tetrahedral mesh. In summary, we have shown that, first, the use of hybrid meshes offers a practical interest as it allows us to accurately solve the neighborhood of the scattering objects while reducing the number of degrees of freedom in the rest of the domain, thus saving computational resources without sacrificing numerical precision. Second, the use of a non-conforming multi-element mesh does yield a notable degradation of the parallel performances as compared to a fully tetrahedral mesh methodology. However, the addition of a constraint at the partitioner level allows us to limit this degradation and to take advantage of the hybrid gains on parallel runs.

Iterative approach with advanced weighting

In the previous section, we have shown that the modification of node and edge weights can lead to improved partitions and better parallel performances. Of course, it seems logical to attach to each mesh cell a metric of computational cost, which depends on the considered physical setting.

Having the perfect weight of each cell, which leads to an almost perfectly balanced partition, is not realistic. This problem can be seen as an optimization problem, with as many parameters as the mesh cells.

In this section, we focus on reducing the computational load imbalance of a DIOGENeS-DGTD simulation by modifying the cell's weight. To do so, we present an iterative approach that can enhance the repartition of weights. At first, all cells and all faces are receiving a weight. This weight can be the same everywhere, or it is possible to set an a priori weight value depending on the simulation settings, as seen in [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method[END_REF]. From here, the mesh is transcripted to its dual graph, as explained in section 6.1.1, and partitioned using MeTiS. Then DIOGENeS-DGTD is run for a few iterations in order to get information concerning the time loop. Local timers are set to recover the computational time of each MPI process. From those times, we can calculate an average simulation time, as well as the deviation from the mean of each process. This allows us to highlight which processes were the most expensive and, as a result, in which subdomain the costly cells are located. This deviation is the key ingredient in the optimization process here. Each node of the graph, representing each element T i , is provided with a weight noted w i . In the partitioning step, the weights are taken into account such that the total weight of the different subdomains are as close as possible. We propose the following weight, which is compatible with any initial weight:

w n+1 i = w n i + m n s , (6.1)
where m n s ∈ R is a metric representing the time deviation in the subdomain s at the iteration n. The methodology is illustrated in figures 6.7 and 6.8. We have applied this weighting strategy to the DIOGENeS-DGTD simulation for the solar cell with the back and front texturing, using the hybrid mesh M2 and running the solver for 100

Strong scalability assessment

We have also performed a strong scalability analysis of the realistic thin-film solar cell presented in chapter 4 by applying the proposed DGTD solver with a tetrahedral mesh of the full solar cell model. This mesh is composed of 1,689,764 elements. This performance analysis is conducted on the Occigen PRACE system hosted by CINES in Montpellier. Each node of this system consists of two Intel Haswell E5-2690@2.6 GHz CPU, each with 12 cores. The parallel speed-up is evaluated for 1000 time iterations of the DGTD-P k solver using a fourth-order low-storage Runge-Kutta time scheme. Performance results are presented in table 6.1 and 6.2. where "Elapsed" is the elapsed time, which is used for the evaluation of the parallel speed-up relatively to the first figure given for each configuration of the DGTD-Pk method. In the first table, the timings include the calculation of an important observable quantity for the considered problem, which is the volumic absorption. The computation of this quantity requires sweep over frequencies in a target spectrum and compute on the fly during the time evolution a discrete Fourier transform of the electrical field, i.e. for each degree of freedom of the DG approximation of the electrical field. The later computation can be performed in a fully parallel way for each element of the mesh. However, the evaluation of the volume absorption is limited to a subvolume, i.e. a layer of the multilayer solar cell model; in the present case, the a-Si:H layer is selected. In figure 6.1, we have used mesh partitions with a uniform weight attached to graph nodes. Thus, we do not take into account the computational load balancing issues raised by this localization of the computation of the volume absorption. Then, the obtained parallel speed-up in table 6.1 is suboptimal and degrades when the interpolation degree k is increased. In table 6.3, we used the iterative partitioning strategy with 20 iterations. The parallel speed-up performed way better, with an improvement from 3.30 to 3.90 for the P 2 interpolation. In addition, the overall computation time is 30% reduced. On the contrary, when the evaluation of the volume absorption is deactivated, the parallel performances are quasi-optimal, i.e.

Hybrid parallelization

Nowadays, modern computers and clusters are based on multi-core CPUs and/or on multi CPU architecture. Such computing node are not totally suited to message passing, as they work on a shared memory. Using MPI to distribute a problem on those cores adds an overhead, mainly due to inter core communications. Such architecture is way more suited to a shared-memory parallelism, which is generally implemented using the OpenMP standard. OpenMP is limited to shared memory systems. However, it offers some advantages. The main one is the incremental programming approach. Indeed, a first parallel version can be obtained by using a few OpenMP directives in the most expensive kernels of the solver. Then, this preliminary parallelization can be incrementally extended to the remaining parts of the solver. Secondly, no communication is done at a multi-core level. Combining OpenMP with MPI can remove both the non-optimal communication on multi-core nodes, and reduce the number of mesh partitions, which naturally mitigate the risk of load imbalance.

In this section, we present preliminary results of such a hybrid parallelization, exploiting both MPI and OpenMP. Our shared memory parallelization within a multi-core node exploits OpenMP tasking, and we have decided to focus on the use of the TASKLOOP construct. OpenMP tasking constructs have been inserted in the most time consuming numerical kernels of the DIOGENeS-DGTD solver. These numerical kernels are concerned with:

• The computation of the right-hand side (RHS) of the system of ordinary differential equations resulting from the DG discretization in the space of system (equation (2.14));

• The update of the physical fields through time-stepping with a Low Storage Runge-Kutta (LSRK) scheme.

The corresponding module components, i.e. subroutines or functions, implement loops over mesh cells with a lot of local matrix-vector products.

Besides the numerical kernels for the computation of the RHS terms and the update of the physical fields through time-stepping with a LSRK scheme, the central computing loop of the DIOGENeS-DGTD solver also involves several auxiliary kernels that may account for a nonnegligible part of the overall computing time of a simulation. Those kernels are concerned with the on-the-fly computation of physical observables and other figures of merit that are required for assessing the results of the simulation. They do not influence the numerical solution, but they cannot be considered as post-processing operations, as most of them need to be updated through the whole simulation. For example, although the DGTD solver computes time-varying distributions of the EM field, most of the physical observables are based on frequency-domain data. Therefore, one needs to compute a discrete Fourier transform of the EM field and related quantities. More importantly, these physical observables are most often associated with particular subvolumes or internal surfaces of the computational domain Ω. Therefore, their computation inherently introduces load balancing issues. All these auxiliary kernels have also been adapted to allow for OpenMP task parallelization.

Performance results

For validating and assessing the performances of the OpenMP-enabled version of the DIOGENeS-DGTD solver, we have selected two test problems of increasing size and complexity: a basic vacuum filled cubic cavity, and a thin-film solar cell. In the following tables, "DGTD-P k " stands for the DGTD solver relying on a uniform polynomial order of degree k. All the numerical experiments reported here have been performed on a single compute node of MareNostrum44 , which is equipped with two sockets Intel Xeon Platinum 8160 CPU with 24 cores, each is running at 2.1 GHz. As a general rule, the number of OpenMP tasks is set to be equal to the number of threads, i.e., of the number of cores since we provide timing measures for runs using between 1 to 32 cores. Our first goal is to assess the achievable scalability with our task-based parallelization of DIOGENeS-DGTD using the OpenMP TASKLOOP construct.

Vacuum-filled cubic cavity.

We make use of a tetrahedral mesh with 20,500 cells. The simulated physical time is set to 10 -8 sec. We recall that for a given mesh when the polynomial interpolation order in the DG formulation is increased, the allowable time-step size, which is subjected to a CFL-like stability condition, decreases. This is why in table 6.4, the number of time-steps increases with the interpolation order. For this model test problem, we obtain speed-up figures that are almost linear up to 16 tasks. As can be expected from the socket architecture of a compute node, the speed-up is lower when we use 32 cores. We make use of a tetrahedral mesh with 41,387 cells. The simulated physical time is set to 10 -15 sec. For this test problem, the computation of the volumic absorption in two of the six constituting material layers of the structure is impacting the achievable speed-up negatively, although the obtained figures are almost perfect up to 8 tasks, as seen in 6.5.

Method

Conclusion

In this chapter, we have presented some profiling of the solver used in this thesis, and some improvements techniques concerning the partitions of the mesh. First, we have shown some imbalances caused by the highly longitudinal aspects of the considered mesh. The extensively expensive surfaces are usually regrouped in the same subdomains due to the common mesh partitioning strategy. This can be improved using weighted graphs, however, this results in a longer global time because of the more expensive communication. The case of hybrid meshes was also considered as it is a complexification, specifically at the non-conforming interfaces. Very bad performances were shown on parallel run with usual partitions. We showed that by avoiding a hybrid interface, splitting results in an almost perfect speed-up. We also gave some information about the simulation of the realistic PV case from chapter 4, and in particular, times of each simulation. In the last part, we showed some results using a hybrid parallelization. This allowed to reduce the stress of the mesh due to much partition. Some tests need to be done for more complex cases, with hybrid meshes and p adaptivity. This work was granted access to the HPC resources of CINES under the allocation 2019-A0060610263 attributed by GENCI (Grand Equipement National de Calcul Intensif) Most of the cases considered in nanophotonics and photovoltaics involve highly heterogeneous materials, with very fine geometric details compared to the size of the domain of interest. Such real-life configurations represent a big challenge for numerical methods because of computational simulations. In a broad sense, discretization methods have to use very fine meshes to represent the smallest geometric details. Moreover, very localized phenomena can lead to high-frequency field distribution which requires the use of refined meshes with also high-order approximation.

In recent years, discretization methods suitable for multiscale problems have been studied extensively. As we have shown in chapter 2 and 3, DG methods are flexible, both from the meshing and from the high-order approximation point of view. It makes it possible to take into account strong heterogeneity of the physical parameters (classically the permittivity ε), as well as very small and curved geometries [START_REF] Viquerat | Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method[END_REF]. From the same family of methods, there is the Hybrid Discontinuous Galerkin (HDG) method presented in [CDL] for Maxwell. One of the key specificity of the HDG method is the splitting of the problem into a global-local formulation. The reason to do so is to make the method globally less expensive by decoupling the volumic degrees of freedom. The counterpart to this is that it requires to solve a linear system at each time-step. Other methods built on top of homogenization techniques such as HMM [START_REF] Henning | A new heterogeneous multiscale method for time-harmonic Maxwell's equations[END_REF] can solve complex heterogeneous problems on coarse meshes. We can also refer to the HHO method [START_REF] Di Pietro | An introduction to Hybrid High-Order methods[END_REF], which shares much similitude in the construction phase with the MHM method studied in the following.

In this chapter, we focus on the MHM method, first presented in [START_REF] Araya | Multiscale Hybrid-Mixed Method[END_REF] for an elliptic model and extented for Darcy's equation in [START_REF] Harder | A hybrid-mixed method for elasticity[END_REF]. It was recently proposed for Maxwell equations in time-domain in [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF]. The fundamental idea of the MHM is to find a way to recover information lost through a rough meshing of the domain. The key ingredients to do so are: (i) the splitting of the initial problem into a global-local formulation; (ii) the construction of the socalled multiscale basis functions that upscale the information from a fine-scale to a coarse scale. Those ingredients together make the MHM method effective, easily parallelizable, and robust since everything is established on a coarse partition.

The efficiency of explicit methods is always suffering when the geometry of the structures requires a high mesh refinement. One possibility is to use local time-step (LTS) methods to split our problems in different regions, each of them using its ideal time-step. The MHM naturally handles this decomposition step. To benefit from it, we develop a new formulation of the MHM method for Maxwell's equation that incorporates the possibility of multiscale time-step (MTS). We extend the results presented in [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF]. In the first part, we present the model and the formulation with a MTS scheme. Then we establish and study energy preservation. Finally, we made two-dimensional numerical tests with comparisons to a standard DGTD method.

Model and preliminaries

In this chapter, we are still considering the system of Maxwell's equations in 3D introduced in chapter 1. We will now start by defining the notations, some of which are specific to the MHM method.

Spatial discretization

The MHM method starts from a discretization of the given spatial domain Ω. We denote by T H the partition of Ω, formed by elements K of maximal size H > 0. From now on, we will denote by T H the coarse mesh or first level mesh and the elements K are denoted by macro element. This subtlety will be important in the second part of this chapter. We also note by ∂T H the skeleton of T H . One can find an illustration of those definitions in figure 7.1. Remark : This partition does not need to fit the geometric details present in Ω.

Time discretization

The second step consists of the coarse time discretization of the time interval [0, T ]. For that we define T ∆t the partition made of the N intervals I n with I n := [t n-1 , t n ] ∈ T ∆t , where N ∈ N * and n ∈ {1, ..., N }. We also define the coarse time-step associated to the interval I n as ∆t n = t nt n-1 > 0.

Remark : This partition may be non-uniform, that is why we keep the exponent n on the time-step ∆t n .

One can find an illustration of those definitions in figure 7.2 We denote by V the broken H 1 space defined as

K 1 K 2 H ∆t H ∆t time space
V := v ∈ L 2 (Ω) : v| K ∈ H 1 (K) for all K ∈ T H (7.1)
with the following scalar product

(u, v) v := K∈T H (u, v) K + K∈T H (∇u, ∇v) K (7.2)
From here we can write the classical system of Maxwell's equations in a local weak formulation as follows.

For each n ∈ {1, ..., N }, find

(E n , H n ) ∈ C 1 (I n , V) × C 1 (I n , V) such that ∀t ∈ I n :        (ε∂ t E n , v) K -(∇ × H n , v) K = (J , v) K , ∀v ∈ V (µ∂ t H n , w) K + (∇ × E n , w) K = 0, ∀w ∈ V E n (t n-1 ) = E n-1 (t n-1 ) H n (t n-1 ) = H n-1 (t n-1 ) (7.3)
with the initial conditions (E 0 (0), H 0 (0)) = (E 0 , H 0 ). We refer to [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF] for the motivation of the use of this space V.

Remark : Please note that ∀n ∈ {1, ..., N }, (E n , H n ) defined in (7.3) are smooth functions of time.

Hybridization and global-local formulation

Next, the continuity of the variables (E, H) is relaxed on the set of boundaries ∂K for all K ∈ T H . Instead, a weak continuity of the tangential component of the electric field E is imposed on ∂T H . We propose the following hybrid weak formulation, originally proposed by Raviart and Thomas in [START_REF] Raviart | Primal Hybrid Finite Element Methods for 2nd Order Elliptic Equations[END_REF] and already used in [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF].

Here we define Λ as the space of the restriction of the tangential component of functions in H(curl; Ω) to the boundaries ∂K i.e.

Λ := v × n K | ∂K ∈ H -1 2 (∂K) for all K ∈ T H : v ∈ H(curl; Ω) (7.4)
with n K the outward normal vector to K. For each

I n ∈ T ∆t , n ∈ {1, ..., N }, find (E n , H n , λ n ) ∈ C 1 (I n ; V) × C 1 (I n ; V) × C 0 (I n ; Λ) such that ∀t ∈ I n :    (ε∂ t E n , v) K -(H n , ∇ × v) K = (J , v) K -(λ n , v) ∂K , ∀v ∈ V (µ∂ t H n , w) K + (∇ × E n , w) K = 0, ∀w ∈ V (ν , E n ) ∂T H = 0, ∀ν ∈ Λ (7.5) with the strong continuity in time E n (t n-1 ) = E n-1 (t n-1
) and H n (t n-1 ) = H n-1 (t n-1 ). We call λ n the hybrid variable on I n . We see that the first two equations of (7.5) are defined locally on each macro element K while the third one is defined on the entire mesh skeleton ∂T H .

Remark : This separation into the global part, and the local part is one of the key aspects of the MHM method. The computational benefits of the method result from this step.

Splitting of electromagnetics fields

We use the linearity of Maxwell's operator to decompose the electromagnetic field (E, H) in two parts: one coming from the external current J, and the other one from the hybrid variable λ. We thus write for all n ∈ {1, ..., N }

E n = E n,J + E n,λ and H n = H n,J + H n,λ (7.6)
where, (E n,J , H n,J ) and (E n,λ , H n,λ ) are respectively solutions of the following local Maxwell problems ∀t ∈ I n :

           ε∂ t E n,J , v K -H n,J , ∇ × v K = (J , v) K , ∀v ∈ V µ∂ t H n,J , w K + ∇ × E n,J , w K = 0, ∀w ∈ V E n,J (t n-1 ) = E n-1 (t n-1 ) H n,J (t n-1 ) = H n-1 (t n-1 ) (7.7) with E n,J (t n-1 ) = E n-1 (t n-1 ) and H n,J (t n-1 ) = H n-1 (t n-1 ) and            ε∂ t E n,λ , v K -H n,λ , ∇ × v K = -(λ n , v) ∂K , ∀v ∈ V µ∂ t H n,λ , w K + ∇ × E n,λ , w K = 0, ∀w ∈ V E n,λ (t n-1 ) = 0 H n,λ (t n-1 ) = 0 (7.8)
With respectively the following initial conditions:

E 0,J (0) = E 0 and H 0,J (0) = H 0 E 0,λ (0) = 0 and H 0,λ (0) = 0 . (7.9)
Remark : We note that the choice of putting all the contributions from the previous I n in the J part as seen in equation (7.9) will have a significant impact on the computational efficiency of the method as it allows to compute problem (7.8) once for all.

Owing to this splitting, the third equation becomes a one-field face formulation for the Lagrange multiplier λ as follows: For each I n ∈ T ∆t , n ∈ {1, ..., N }, find λ n ∈ Λ such that ∀t ∈ I n :

ν , E n,λ ∂T H = -ν , E n,J ∂T H , ∀ν ∈ Λ, (7.10)
We precise that when λ is used in the exponent of E n,λ , one should see E n,λ n . The additional exponent is omitted here for the sake of readability.

One-level MHM discretization

We start with the space-time discretization of the Lagrange multiplier λ. We called this step the one-level or first-level discretization.

Discretization of the hybrid variable

We choose a finite dimensional space to approximate the variables in C 0 (I n ; Λ). We are especially looking for λ H (t) in such a space satisfying (7.10). In this work, we choose to take the Lagrange multipliers λ H constant in time per time-slab I n i.e. λ H ∈ P 0 (I n ), n ∈ {1, ..., N } and polynomial on T H , i.e. λ H ∈ P 0 (I n ) ⊗ Λ H with:

Λ H := v H ∈ Λ : v H | F ∈ P m (F ) 3 for all F ∈ ∂T H , m ∈ N, (7.11)
and look for λ n H in Λ H . Similarly, equation (7.6) becomes a discrete decomposition as:

E n H = E n,J + E n,λ H and H n H = H n,J + H n,λ H (7.12)
From here equation (7.8) becomes:

Find (E n,λ H , H n,λ H , λ n H ), n ∈ {1, ..., N } in (V, V, Λ H ) such that ∀t ∈ I n :            ε∂ t E n,λ H , v K -H n,λ H , ∇ × v K = -(λ n H , v) ∂K , ∀v ∈ V µ∂ t H n,λ H , w K + ∇ × E n,λ H , w K = 0, ∀w ∈ V E n,λ H (t n-1 ) = 0 H n,λ H (t n-1 ) = 0 (7.13)
with the initial conditions: E 0,λ H = 0 and H 0,λ H = 0 (7.14)

and the global problem: For each ), n ∈ {1, ..., N }, for all K × I n ∈ T H × T ∆t , as the solution of the following problem:

I n ∈ T ∆t , n ∈ {1, ..., N }, find λ n ∈ Λ such that ∀t ∈ I n : ν , E n,λ H ∂T H = -ν , E n,J ∂T H , ∀ν ∈ Λ H , (7.15
           ε∂ t η n,E i , v K -η n,H i , ∇ × v K = -ψ i n • n K , v ∂K , ∀v ∈ V(K) µ∂ t η n,H i , w K + ∇ × η n,E i , w K = 0, ∀w ∈ V(K) η n,E i (t n-1 ) = 0 η n,H i (t n-1 )
= 0 (7.17) with V(K) the space of functions in V restricted to K. These functions are called the multiscale basis function. They serve as a basis to reconstruct E n,λ H and H n,λ H . Multiplying each equation by β n i , we obtain:

           εβ n i ∂ t η n,E i , v K -β n i η n,H i , ∇ × v K = -β n i ψ i n • n K , v ∂K , ∀v ∈ V(K) µβ n i ∂ t η n,H i , w K + ∇ × β n i η n,E i , w K = 0, ∀w ∈ V(K) η n,E i (t n-1 ) = 0 η n,H i (t n-1 ) = 0 (7.18)
and we sum for each i in {1, ..., dim Λ H }:

                     ε dimΛ H i=1 β n i ∂ t η n,E i , v K - dimΛ H i=1 β n i η n,H i , ∇ × v K = - dimΛ H i=1 β n i ψ i n • n K , v ∂K , ∀v ∈ V(K) µ dimΛ H i=1 β n i ∂ t η n,H i , w K + ∇ × dimΛ H i=1 β n i η n,E i , w K = 0, ∀w ∈ V(K) η n,E i (t n-1 ) = 0 η n,H i (t n-1 )
= 0 (7.19) From here and using equations (7.16) and (7.13) we can identify ∀t ∈ I n :

E n,λ H = dimΛ H i=1 β n i η n,E i and H n,λ H = dimΛ H i=1 β n i η n,H i , (7.20)
with the initial conditions η n,E i (t n-1 ) = 0 and η n,H i (t n-1 ) = 0. Using equation (7.20) in (7.12), we obtain ∀t ∈ I n :

E n H = E n,J + dimΛ H i=1 β n i η n,E i and H n H = H n,J + dimΛ H i=1 β n i η n,H i (7.21)
From there we note the set of (η n,E i , η n,H i ) with i ∈ {1, ..., dim Λ H }, the multiscale basis functions. They are upscaling local information from the problems (7.17) to the coarse mesh. One can see an example in figure 7.3. Let us point out that the set of β n i is a fundamental part of the resulting algorithm. Indeed, (β n i ) are the only globally coupled d.o.f. in the MHM method and are solution of a linear system defined on the coarse mesh skeleton ∂T H . This linear system results from (7.10) and decomposition (7.20) as for each n ∈ {1, ..., N } find β n i such that ∀t ∈ I n :

ν H , dimΛ H i=1 β n i η n,E i ∂T H = -ν H , E n,J ∂T H , ∀ν H ∈ Λ H , (7.22)
given η n,E i and E n,J .

η n,E i d i If one assumes that close analytical formula are available for all n ∈ {1, ..., N } for (η n,E , η n,H ) and (E n,J , H n,J ), it is still necessary to obtain the d.o.f. β n i to compute the solution from system (7.22). However, for real application test case, we do not have access to analytical formulas, meaning a second level of discretization is necessary to make the method practical. This is the subject of the next section.

Two-level MHM method

As stated at the end of the first section, analytical formulas for both η n,ξ and ξ n,J , ξ ∈ {E, H} are not available for realistic cases. In order to solve the second level problems defined in (7.7) and (7.13), we choose to use a spatially centered DG scheme coupled with a second order Leap-Frog (LF) time scheme.

DG on each macro element K

The second level of space discretization will be called local mesh or fine mesh, unlike the global mesh or coarse mesh T H . More precisely, in each macro element K of T H we build a fine partition noted T K h , formed by elements κ of characteristic length h K > 0. We precise that when h is used in the indices of T K h , one should see T K h K . The additional exponent is omitted here for the sake of readability. One can find an illustration of this definition in figure 7.4. We also define the set of faces of the local mesh T K h by F K h as well as the set of faces on the boundary of K as F K 0 .

Remark : It is important to notice that no continuity at the junction between macro elements is mandatory, as illustrated in figure 7. 4.

T H Ω T K 2 h K 2 T K 1 h K 1 Figure 7.4 | Space discretization. Second level space discretization. Each macroelement K are discretized in submeshes T K h K .
We discretize local problems equations (7.7) and (7.17) with respect to the space variable by selecting a local finite dimensional space V h (K) whose functions are defined over {T K h } h>0 . The two discrete problems read as follows:

1. Find for n ∈ {1, ..., N }, (E n,J h , H n,J h ) in V h (K) × V h (K) such that ∀t ∈ I n :            ε∂ t E n,J h , v h K -H n,J h , ∇ × v h K = (J , v h ) K , ∀v h ∈ V h (K) µ∂ t H n,J h , w h K + ∇ × E n,J h , w h K = 0, ∀w h ∈ V h (K) E n,J h (t n-1 ) = E n-1 h (t n-1 ) H n,J h (t n-1 ) = H n-1 h (t n-1 ), (7.23)
with E 0,J h = E 0 and H 0,J h = H 0 .

Find for

n ∈ {1, ..., N }, (η n,E i,h , η n,H i,h ) in V h (K) × V h (K), such that for all i ∈ {1, ..., dimΛ H } such that ∀t ∈ I n :              ε∂ t η n,E i,h , v h K -η n,H i,h , ∇ × v h K = -ψ i n • n K , v h ∂K , ∀v h ∈ V h (K) µ∂ t η n,H i,h , w h K + ∇ × η n,E i,h , w h K = 0, ∀w h ∈ V h (K) η n,E i,h (t n-1 ) = 0 η n,H i,h (t n-1 ) = 0. ( 7
ν h , dimΛ H i=1 β n i η n,E i,h ∂T H = -ν h , E n,J h ∂T H , ∀ν h ∈ Λ H , (7.25)
From here we can reconstruct our fields:

E n H,h = E n,J h + dimΛ H i=1 β n i η n,E i,h and H n H,h = H n,J h + dimΛ H i=1 β n i η n,H i,h (7
.26)

LF2 time integration on each time-slab I n

In this part, we devise the fully discrete two-level MHM method by setting up a time-marching scheme for the local problems equations (7.23) and (7.24). A second-order Leap-Frog scheme with a multiscale time-stepping (MTS) strategy is proposed.

The second level time discretization comes from a decomposition of each I n in a uniform partition of time-step ∆τ n,K > 0 on each K ∈ T H , i.e. t n = t n,0 < t n,1 < .... < t n,M n,K = t n+1 , with M n,K ∈ N, where t n,k = t n,0 + k∆τ n,K , k ∈ {1, ..., M n,K }. In the same way we call micro time-slab or fine time-slab the interval I n,m = [t n,m , t n,m+1 ] and ∆τ n,K the micro time-step associated to I n,m .

Remark : Note that this partition can be different in each element K ∈ T h . Strictly speaking, we should have one more exponent related to the macro element since the fine discretizations are independent from one K to another.

Next, we focus on the intrinsic MTS of the MHM method within a simple framework by assuming that T ∆t is uniform with characteristic length ∆t n = ∆t, n ∈ {1, ..., N }. We choose to divide each I n ∈ T ∆t of each K into M K intervals t n = t n,0 < t n,1 < . . . < t n,M K -1 < t n,M K = t n+1 , so that we can omit the time exponent from ∆τ n,K i.e. ∆τ n,K = ∆τ K , n ∈ {1, ..., N }. On K, we denote by E n,m H,h (resp. H n,m H,h ) the field that approximate E n H,h (resp. H n H,h ) at the discrete local time t n,m , m ∈ {0, ..., M K }. Since t n,M K = t n,0 + M K ∆τ K = t n,0 + ∆t n = t n+1,0 , we keep the continuity of E n,M H,h = E n+1,0 H,h and H n,M H,h = H n+1,0 H,h between each macro time-slab as illustrated in figure 7.5. Owing to these definitions, the centered DG method applied to equation (7.24) reads as follows.

K 1 K 2 H ∆t ∆τ K 2 H ∆t ∆τ K 1 time space
Let K ∈ T H , find for each m ∈ {1, . . . , M K }, (E n,m,λ H h , H n,m,λ H h ) ∈ C 1 (I n , V h (K))×C 1 (I n , V h (K)), for all I n ∈ T ∆t , such that                      1 ∆τ K εη m+ 1 2 ,E i,h , φ K k K = 1 ∆τ K (εη m-1 2 ,E i,h , φ K k K + η m,H i,h , ∇ × φ K k K -( η m,H i,h , φ K k ) F K 0 -(ψ i , φ K k ) ∂K , 1 ∆τ K µη m+1,H i,h , φ K k K = 1 ∆τ K (µη m,H i,h , φ K k K -η m+ 1 2 ,E i,h , ∇ × φ K k K + ( η m+ 1 2 ,E i,h , φ K k ) F K h , (7.27) with η 1 2 ,E i,h = 0 and η 0,H i,h = 0, for all (v h , w h ) ∈ V h (K) × V h (K).
One can retrieve the definitions of the DG operators • and • in chapter 2.

When applied to (7.7), it becomes: For each

K ∈ T H , m ∈ {0, . . . , M K } find (E n,m,J h , H n,m,J h ) ∈ V h (K) × V h (K), for all I n ∈ T ∆t , such that          ε E n,m+ 1 2 ,J h -E n,m-1 2 ,J h ∆τ K , v h T K h -H n,m,J h , ∇ × v h T K h + ( H n,m,J h , v h ) F K 0 = (J n,m , v h ) T K h , µ H n,m+1,J h -H n,m,J h ∆τ K , w h T K h + E n,m+ 1 2 ,J h , ∇ × w h T K h -( E n,m+ 1 2 ,J h , w h ) F K h = 0 , (7.28) given E n, 1 2 ,J h = E n-1 H,h and H n,0,J h = H n-1 H,h with E 0, 1 2 ,J h = E 0 and H 0,0,J h = H 0 , for all (v h , w h ) ∈ V h (K) × V h (K).
Remark : As one can see, there is no more exponent n in the multiscale basis problem (7.27). This invariance with respect to the macro time-slabs comes from two points: (i) the right-hand side of the problem (7.27), which only depends on the basis function ψ; (ii) the macro timeslab initial conditions are (0, 0). From there, one can solve the problem (7.27) once for all

I n i.e. I n,m K = [t n,m K , t n,m+1 K
]. We choose to not write explicitly this dependence if there is no ambiguity.

With regard to the global problem (7.25), a choice must me made for its resolution. In [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF], where no MTS strategy is envisaged, the global problem is solved once per time-slab I n , n ∈ {1, ..., N }, more precisely at t n,M K . In our case, we envisage the evolution through one macrotime slab. In order to take into account history inside each I n , we propose the following equivalent formulation: For all

I n ∈ T ∆t , find λ n H ∈ Λ H such that ν H , ˆIn dimΛ H i=1 β n i η n,E i,h ∂T H = -ν H , ˆIn E n,J h ∂T H , ∀ν H ∈ Λ H , (7.29)
Remark : For the moment, equation (7.29) is not discretized in time. Different integration quadrature are possible. We show in the following that the choice of the integration quadrature affects the stability of the scheme.

Energy preservation

The global problem is approximated for all n ∈ {1, ..., N } by:

∆τ K M m=1 E n,m+ 1 2 H,h + E n,m-1 2 H,h 2
K) = ε 2 E n,m+ 1 2 H,h (
0,T K h + 1 2 µH n,m H,h , H n,m+1 H,h T K h . ( 2

.31)

The first level energy given on I n at discrete time t n by

E n H,h = K∈T H E n,M K +1/2 H,h (K) (7.32)
One has the following energy , ∀n ∈ {1, ..., N }:

E n H,h = K∈T H E 0,1/2 H,h (K) + K∈T H n k=0 M K m=1 (J k,m , E k,m+ 1 2 + E k,m-1 2 2 ) T H ,
So that, if one has a vanishing current J, the energy is preserved over discrete times.

Let n ∈ {0, .., N }. We write the local problems on each I n × K and use the equivalent discrete equations resulting from (7.27), (7.28) and (7.29). Given (E

n, 1 2 ,λ H H,h , H n,0,λ H H,h ), find for each m ∈ {1, . . . , M }, (E n,m+1/2 H,h , H n,m H,h ) ∈ C 1 (I n , V h (K)) × C 1 (I n , V h (K)) such that              ε E n,m+ 1 2 H,h -E n,m-1 2 H,h ∆τ K , v h T K h -H n,m H,h , ∇ × v h T K h + ( H n,m H,h , v h ) F K 0 = -(λ n H , v h ) ∂K , +(J n,m , v h ) T K h , µ H n,m+1 H,h -H n,m H,h ∆τ K , w h T K h + E n,m+ 1 2 H,h , ∇ × w h T K h -( E n,m+ 1 2 H,h , w h ) F K h = 0 , (7.33) with E n, 1 2 H,h = E n-1,M K + 1 2 H,h , H n,0 H,h = H n-1,M K H,h
and J n,m the external current computed at time t n,m .

Remark : For the second equation, on faces of F K h \ F K 0 i.e. the internal faces, we use

E n,m+ 1 2 H,h = E n,m+ 1 2 H,h and
w h = w h .
Proof. One has for n ∈ {1, ..., N } and m ∈ {1, ..

., M K } E n,m+1/2 H,h (K) -E n,m-1/2 H,h (K) = ε 2 E n,m+ 1 2 H,h 2
0,T K h + 1 2 µH n,m H,h , H n,m+1 H,h T K h - ε 2 E n,m-1 2 H,h 2
0,T K h - 1 2 µH n,m-1 H,h , H n,m H,h T K h = 1 2 ε (E n,m+ 1 2 H,h -E n,m-1 2 H,h ), E n,m+ 1 2 H,h + E n,m-1 2 H,h T K h + 1 2 µ(H n,m+1 H,h -H n,m-1 H,h , H n,m H,h T K h .
Adding the equations (7.28) to its λ counterpart gives us this complete form of system (7.28):

                 ε E n,m+ 1 2 H,h -E n,m-1 2 H,h ∆τ K , v h T K h = H n,m H,h , ∇ × v h T K h -( H n,m H,h , v h ) F K 0 -(λ n H , v h ) ∂K + (J n,m , v h ) T K h , µ H n,m+1 H,h -H n,m H,h ∆τ K , w h T K h = -E n,m+ 1 2 H,h , ∇ × w h T K h + ( E n,m+ 1 2 H,h , w n,m H,h ) F K h , (7.34)
Testing with respectively the fields

E n,m+ 1 2 H,h +E n,m-1 2 H,h 2 and H n,m
H,h , and the second equation of (7.34) at discrete step m -1, with H n,m H,h , one obtains the three equations:

                                                               ε E n,m+ 1 2 H,h -E n,m-1 2 H,h ∆τ K , E n,m+ 1 2 + E n,m-1 2 2 T K h = H n,m H,h , ∇ × E n,m+ 1 2 + E n,m-1 2 2 T K h -( H n,m H,h , E n,m+ 1 2 + E n,m-1 2 2 ) F K 0 -(λ n H , E n,m+ 1 2 + E n,m-1 2 2 ) ∂K +(J n,m , E n,m+ 1 2 + E n,m-1 2 2 ) T K h , µ H n,m+1 H,h -H n,m H,h ∆τ K , H n,m H,h T K h = -E n,m+ 1 2 H,h , ∇ × H n,m H,h T K h +( E n,m+ 1 2 H,h , H n,m H,h ) F K h , µ H n,m H,h -H n,m-1 H,h ∆τ K , H n,m H,h T K h = -E n,m-1 2 H,h , ∇ × H n,m H,h T K h +( E n,m-1 2 H,h , H n,m H,h ) F K h , (7.35) This gives,                                                ε E n,m+ 1 2 H,h -E n,m-1 2 H,h ∆τ K , E n,m+ 1 2 + E n,m-1 2 2 T K h = H n,m H,h , ∇ × E n,m+ 1 2 + E n,m-1 2 2 T K h -( H n,m H,h , E n,m+ 1 2 + E n,m-1 2 2 ) F K 0 -(λ n H , E n,m+ 1 2 + E n,m-1 2 2 ) ∂K +(J n,m , E n,m+ 1 2 + E n,m-1 2 2 ) T K h , µ H n,m+1 H,h -H n,m-1 H,h ∆τ K , H n,m H,h T K h = -E n,m+ 1 2 H,h + E n,m-1 2 H,h , ∇ × H n,m H,h T K h +( E n,m+ 1 2 H,h + E n,m-1 2 H,h , H n,m H,h ) F K h , (7.36
) Thus, one can estimate the variation of the fully discrete energy over I n,m as:

E n,m+1/2 H,h (K) -E n,m-1/2 H,h (K) ∆τ K = H n,m H,h , ∇ × E n,m+ 1 2 + E n,m-1 2 2 T K h -( H n,m H,h , E n,m+ 1 2 + E n,m-1 2 2 ) F K 0 -(λ n H , E n,m+ 1 2 + E n,m-1 2 2 ) ∂K +(J n,m , E n,m+ 1 2 + E n,m-1 2 2 ) T K h , -E n,m+ 1 2 H,h + E n,m-1 2 H,h , ∇ × H n,m H,h T K h +( E n,m+ 1 2 H,h + E n,m-1 2 H,h , H n,m H,h ) F K h (7.37)
Summing (7.37) for m ∈ {1, ..., M } make a lot of terms telescoping , which at the end gives us:

E n H,h -E n-1 H,h = K∈T H M K m=1 ∆τ K (E n,m+1/2 H,h (K) -E n,m-1/2 H,h (
K)) = K∈T H M K m=1 -∆τ K (λ n H , E n,m+ 1 2 + E n,m-1 2 2 ) ∂K + K∈T H M K m=1 ∆τ K (J n,m , E n,m+ 1 2 + E n,m-1 2 2 ) T K H , (7.38)
Finally, if we write the global evolution of the discrete energy over [0, t n ], we obtain,

E n H,h -E 0 H,h = n k=0 E k+1 H,h -E k H,h = n k=0 K∈T H M K m=1 -∆τ K (λ n H , E k,m+ 1 2 + E k,m-1 2 2 ) ∂K + n k=0 K∈T H M K m=1 ∆τ K (J k,m , E k,m+ 1 2 + E k,m-1 2 2 ) T K H , (7.39)
where we have used the continuity of the field on the boundary of I n . We now use the global problem (7.30) on each I n with test function µ H = λ n H . We find:

E n+1 H,h -E 0 H,h = K∈T H n k=0 M K m=1 ∆τ K (J k,m , E k,m+ 1 2 + E k,m-1 2 2 ) T K H ,

Stability

In order to prove that our scheme is stable, we study if

E n,M +1/2 H,h
is a positive definite quadratic form of E n,m+ 1 2 and H n,m . We suppose that each T K h is a quasi uniform mesh, i.e. T K h is shape regular and ∃ζ We recall here inverse inequalities [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Ern | Evaluation of the condition number in linear systems arising in finite element approximations[END_REF]. We can assume that there exists a constant C > 0 such that for all p ∈ P l (T K h ), H,h 0,τ j H n,m H,h 0,τ i , moreover, using quasi-uniformity, one as: As a result, we have shown that the energy is a bounded positive quadratic form which proves the stability.

K > 0 ∀τ i ∈ T K h , h τ i ≥ ζ K h K with h K = max τ i ∈T K h h K τ i .
∇ × p 2 0,τ ≤ Ch -1 τ p 2 0,τ ∇ × p 2 0,T K h = τ ∈T K h ∇ × p 2 0,τ ≤ τ ∈T K h Ch -1 τ p 2
∆τ K E n,m+ 1 2 H,h , ∇ × H n,m H,h T K h ≤ ∆τ K C 1 ζ K h K τ ∈T K h E n,m+ 1 2 H,h 0,τ H n,m H,h 0,τ |∆τ K E n,m+ 1 2 H,h , H n,m H,h T K h | ≤ ∆τ K C 2 C 3 2ζ K h K τ i ∈T K h τ j ∈V i E n,m+ 1 2 H,h 0,τ i H n,m H,h 0,τ i + ∆τ K C 3 C 4 2ζ K h K τ i ∈T K h τ j ∈V i E n,
- ∆τ K C 1 ζ K h K τ ∈T K h E n,m+ 1 2 H,h 0,τ H n,m H,h 0,τ - ∆τ K C 2 C 3 2ζ K h K τ i ∈T K h τ j ∈V i E n,m+ 1 2 H,h 0,τ i H n,m H,h 0,τ i - ∆τ K C 3 C 4 2ζ K h K τ i ∈T K h τ j ∈V i E n,
- ∆τ K C 2 C 3 C 5 2ζ K h K τ i ∈T K h τ j ∈V i ε 2 E n,m+ 1 2 H,h 2
0,τ i + µ 2 H n,m H,h 2 0,τ i - ∆τ K C 3 C 4 C 5 2ζ K h K τ i ∈T K h τ j ∈V i ε 2 E n,m+ 1 2 H,h 2
0,τ j + µ 2 H n,m H,h 2 0,τ i = 1 2 - ∆τ K C1 2ζ K h K - ∆τ K C2 4ζ K h K - ∆τ K C3 4ζ K h K τ ∈T K h ε E n,

Numerical results

This section presents preliminary results obtained with a MHM-DGTD code with MTS strategy. The Fortran90 code solves the Maxwell's equations in 2D.

Algorithm

First of all, even if the results presented in this part are obtained using a non parallel implementation of the MHM-DGTD method, we want to highlight all the potential of the method in the following description of the algorithm steps of algorithm 3 (see below).

Algorithm 3 Summary of the MHM method with MTS algorithm 1: Set E 0 K and H 0 K , ∀K ∈ T H from initial conditions 2: for n = 1 → N do Set E n,J K = E n-1 K and H n,J K = H n-1 K , ∀K ∈ T H 3:

for k ∈ T H do Loop on macro elements Solve the global problem and update E n,0 K and H n,0 K , ∀K ∈ T H 15: end for Compared to DG methods for which all the communications between cells are done via the numerical fluxes, the MHM method is almost entirely local if we neglect the cheap global problem in line 14. Each current problem defined in equation (7.28) is associated with a macro element K. They are all considered as single independent problems with their own initial and boundary conditions changing over the time-slabs. As a result, we expect an almost perfect scalability with a parallelization of the local problems. Indeed, tasks 6-7 and 9-10, which represent respectively local problems (7.27) and (7.28), are embarrassingly parallel, i.e. that there is no dependency between tasks.

The aforementioned parallel strategy is the intuitive way of thinking the parallelism of the MHM. However, more advanced possibilities are considerable, as for example, a two-level parallelization or a hybrid parallelization could allow to accelerate computationally expensive local problems.

Since the MHM method involves two levels of discretization, all the standard numerical parameters are doubled. The two mesh sizes H and h, as well as the interpolation order of the basis functions, both for the first level space λ H , denote by p e , and for the second level space V h denoted by p k . As shown by Harder et al in [START_REF] Harder | A hybrid-mixed method for elasticity[END_REF] for elasticity, in the case where H = h, the couple (p e , p k ) must verify p k = p e + 1 if p e is odd and p k = p e + 2 if it is even. We will see that this will influence our results as well.

Convergence of the degenerated scheme

In this subsection, we start with a convergence study of the degenerated time scheme, i.e. when there is no second level in time (∆τ K = ∆t, ∀K ∈ T H ). For this purpose, we consider the analytical solution of the cavity mode given on a domain Ω = [0, 1] 2 with Perfect Electric Conductor (PEC) • Global convergence analysis (figures 7.9)

Here, we look at the classical convergence, i.e. when (H, h K , ∆t, ∆τ K ) → 0 simultaneously. The pair of order used here is still (p e = 1, p k = 3), with a spatial refinement H h going from 1 to 8. All the four parameters are linked to each other from the following relations: i) H h = {1, 2, 4, 8} from subfigures (a) to (d); ii) ∆τ K = min K∈T H CF L h K ; iii) ∆t = α∆τ K with α ∈ {1, 5, 19, 15, 20}, represented by the different color lines. Thus setting only one of the four parameters to 0 amounts setting all of them to 0. As previously, we observe an order 2 for the degenerated case i.e. ∆t = ∆τ . Focusing first on the subfigure (a), we observe an order 1. However, as H decreases we remark a change of order from 2 to 1. This does not seem to happen in subfigure (d) but this might be due to the fact that H h = 8 and we might observe the same behavior for H smaller than the values achieved (here 1.56e -2 ). From this, one would infer min(p e + 1, p K , e, k) = 1.

Remark : It is difficult to conclude if it is a loss of order or a limitation of the scheme as we are expecting at best a second order. • Convergence analysis in terms of the global parameters (∆t, H) (figure 7.10) The third case is motivated by the two levels of the MHM method. One can ask how the first level behaves when the error of the second level is negligible. For this case we study the convergence in (∆t, H) by decoupling the errors using (∆τ K , h K ) small enough so as not to pollute the errors in (∆t, H). We look at different ratios ∆t H in subfigures (a),(b) and (c). We use a couple of orders (p e , p K ) = (1, 2) as we avoid the case h K = H. The different lines represent several values of (∆τ K , h K ) with ∆τ K = CF L h for each of these ratios. In every tests the order obtained is 2 (or even 2.5). We infer that min(e, p e + 1) = 2, so e ≤ 2. In conclusion, we have shown numerically that the MHM method with MTS scheme converges. First of all, we confirmed its unconditional stability numerically with respect to the macro time and space parameters, by showing the energy conservation and its positivity for any (∆t, H). Then we tried to investigate the convergence order of the respective parameter. However, one test shows a decrease from 2 to 1 in the convergence order for specific parameters, especially when the number of local time-steps per global time-step increases. The increasing complexity of local problems can explain this. Indeed, the increase of ∆t amounts to increasing the final time of each local problems. We know that doing so increases the error, especially the dispersion one. Therefore it makes the local problems more difficult to solve. To keep the same accuracy, one should refine the local problems. This is a comportment we saw in figure 7.9.

Since for the moment, we cannot rely on a rigorous error analysis, we can only infer that the error estimate might be more intricate than O(H p e +1 + h p k + ∆τ k + ∆t e ). Indeed, figure 7.9 might show that a complex relationship between the two levels must be respected to reach the optimal convergence, and even maybe even a super-convergence as in figure 7.10. In this part, we aimed at comparing the MHM method and the DGTD method on a realistic multiscale nanophotonic case. To do so, we select a prototype problem of a photonic crystal structure tested in [START_REF] Lanteri | The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media[END_REF]. The nanophotonic device is composed of a silica chip filled with 15 round silicon holes of radius 0.3125. For this study, we place the device in a domain Ω = [0, 10] × [0, 10] made of air as shown in figure 7.11. A monochromatic plane wave of frequency 0.15 (dimensionless frequency) is coming from the left of the structure. The permittivity used for the 3 materials are ε = 1 for the surrounding air, ε = 1.5 for the silica chip and ε = 3.14 for the silicon holes. Absorbing Boundary Condition (ABC) are used on the whole boundary ∂Ω.

Numerical features

To be able to compute a numerical error for this realistic case, we start by computing one reference solution. To do so, we use DGTD method on a highly refined mesh with P 2 basis functions. The reference solution is plot in figure 7.14(a). We choose to compute the error on nine fixed points, visible in the figure 7.12, instead of projecting the solutions on the same mesh. The study is carried out as follows: i) We take a very coarse mesh M1, not fitting to the geometry (see in figure 7.13(a) that M1 is not meshing the spheres). We show that the MHM on this mesh obtains a correct error. ii) More and more refined meshes are used for the DGTD method until the MHM error is more or less reached. We used three meshes as shown in figure 7.13.

Results

The information related to the three mesh sizes used is presented in the table 7.1 as well as the maximum error at the control points. One can see that on the very coarse mesh, the MHM method is performing way better than the DGTD method. This is predictable given the absence of the silica holes in M1. However, to reach a solution close to the MHM results, we use a mesh with nearly 40 times more d.o.f.. Of course, we omit the second level discretization deliberately. However, we can see from the accuracy of the visualization of figure 7.14(b), that the structure is thinly discretized.

Computation times

Now that we have shown that the MHM can precisely solve a complex problem on a very coarse mesh, we will look at the performance of the method in terms of computation times.

For each simulation, we retrieve the times of the different computation phases in both solvers, i.e. the computation of the multiscale basis functions η ξ before the time-loop, the resolution of local problems to retrieve (E J , H J ) at each time-step, the resolution of the global problem, and the summation E = E J + E λ . The two firsts are referred to as local times while the two lasts are referred to as the global times.

In figure 7.15, we present the comportment of the MHM method in regards to the macro timestep used. Fine time-steps are set thanks to the micro CFL, and we increase ∆t such that the ratio ∆t/∆τ K goes from 1 to 8. We observe a quasi-optimal scalability in terms of global times.

Remark : We do not compare computational times between MHM and DGTD. Indeed: (i) the MHM targets in particular potentially unsolvable problems for the DGTD method and should com- pete on those type problems; (ii) the MHM needs to be parallelized in order to compete; and (iii) it is necessary to understand the interplay between the parameters (T H , T ∆t , p e , p K ) in order to use the method as efficiently as possible.

Conclusion

In this chapter we presented a MHM method using a DGTD method at the second level. We proposed a MTS time scheme to benefit from the independence of local problems. We proved the preservation of the energy and the stability of the method using a trapezoidal rule for the computation of the global problem, and we illustrated it numerically. However, a theoretical study of the convergence should help to understand the behavior of the error with MTS scheme. Moreover, the use of higher-order time scheme is also an objective for the future of the explicit MHM.

We also showed the capacity of the MHM method to solve accurately complex, realistic problems on very coarse meshes. However, for the moment, the usefulness of the method is limited. Indeed, its high computational cost requires an adapted parallelism to be competitive or at least reasonable in terms of computation time. Several levels of parallelism can be considered, with the most logical first approach being the parallelization of local problems. One can also think about using parallel DGTD solver for the local problems.

We also think that the method can be well suited for implicit time scheme and for frequencydomain problems [START_REF] Chaumont-Frelet | A multiscale hybrid-mixed method for the Helmholtz equation in heterogeneous domains[END_REF]. Indeed, the bottleneck of these methods for real test cases is, most of the time, the amount of memory available to store the linear system. However, in the same way, that HDG method did, MHM method is spreading out the computation load in the local problems whereas the global linear system to solve involves only the skeleton of the coarse mesh. The purpose of this last chapter is to summarize the content of this manuscript and to identify possible extensions.

Summary

This thesis aimed at proposing a 3D discontinuous Galerkin time-domain method able to simulate complex and realistic PV devices. First, an introduction to the context of photovoltaics, and its link with electromagnetism has been made. A presentation of Maxwell's equations, and its coupling with the generalized dispersive model was done. This dispersion model can represent accurately a wide diversity of dispersive materials on a broad band, and especially semiconductors which are one of the main components in solar cells. Then we have presented the DGTD method with its spatial discretization. Time discretization using RK and LSRK schemes has also been introduced with some words on stability. The imposition of a plane wave in the computational domain has been discussed, as well as the decomposition of the domain in a total field region and a scattered field region. Finally, the domain truncature has been discussed, with the introduction of PMLs, PBCs and the other classical boundary conditions.

In the third chapter, we have presented two methodological contributions that aimed to improve the performance and the accuracy of the DGTD method in the particular for solar cell problems. In a first part, the possibility to use local approximation order was investigated. The selected order repartition is based on the number of points per wavelength, which is required to obtain a correct approximation of wave phenomena. This strategy was applied on a realistic silicium solar cell. Both accuracy of the approximation (< 1% relative error) and sequential speed-up (≈ 2) was interesting. The second contribution consist of using hybrid non-conforming meshes. This strategy allowed the use of a refined unstructured mesh at the material junctions while having a structured coarse grid in the homogeneous regions. This strategy has been tested and validated on simple test cases and on the realistic solar cell. Coupling both strategies showed an even bigger speed-up (≈ 3).

In the fourth chapter, we tried to demonstrate the capabilities of our DGTD solver on realistic photovoltaic cases. In the first part we presented state-of-the-art methods to deal with such applications. The most important quantities of interest were presented (Absorption, EQE, Jsc). Then, we showed how to deal with geometrical measures obtained in laboratory. Different pre-and post-processing were presented (smoothing raw measures, building periodic meshes from random textures) using methods taken from digital image processing. After that, we have presented the simulations of different type of solar cells using DGTD. The first one comes from a collaboration with C2N and IPVF on GaAs based solar cells. We compared the results obtained from different method. The second case dealt with solar cells based on conical gratings. Those particular cells enhanced a lot of resonances, which required a long physical time to catch. The last one consists of a very large solar cell with random textures obtained by AFM. We showed the capability of the DGTD solver to deal with large geometries (10µm 2 ) and big meshes (more than one million cells).

In the fifth chapter, we coupled our DGTD solver with an optimization method called Efficient Global Optimization (EGO). This method is based on a metamodelling, which aims to replicate the objective function. This method is particularly well adapted to optimize expensive and unknown (black-box) objective function. Several optimizations were conducted. The first one is the enhancement of color rendering of a solar cell, for which we were able to provide a nanotexturation that reduces and unifies the parasitic reflections, especially in the red / near infrared region. The second case is the optimization of a pyramid based nanotexturation. The main idea here was to maximize the density of short-circuit current while keeping a fixed quantity of absorber (InGaAs). The q-EGO method was used for this case, which allow parallel evaluation of the objective function per optimization step. The optimized design here increased the short-circuit current density by 23%. Load unbalances were observed during the optimization loop, especially when the configuration present huge differences.

In the sixth chapter, the parallel performances of our code has been assessed. In a first part we presented the parallel implementation of our solver, which relied on a partition of the mesh. A numerical experiment with the Metis mesh partitioning tool and an uniform weighting, speedups and efficiencies of our parallel MPI implementation were assessed on a realistic solar cell case. Bothh hybrid and tetrahedral meshes were tested and compared. Results showed that this simple weighting worked better for tetrahedral meshes. For a good scaling with hybrid meshes, we added a constraint for the partitioning of hybrid faces, which allowed us to solve this issue. We also proposed an iterative partitioning algorithm, which used ongoing parallel simulation to weight the different mesh cells. This algorithm allows to balance the load on each CPU, especially when the considered case has complex features (PMLs, TFSF, observable computations). On a realistic solar cell, the load unbalance went from 15% to 2%. Preliminary results on hybrid MPI/openMP parallelization were presented.

In the last chapter, we have presented of a recent method. The MHM is particularly well adapted to deal with multiscale problems with geometric details on coarse mesh. To be used for realistic applications, the MHM needs a second numerical method to solve local problems in the element of the coarse partition. We chose to use DGTD with centered flux. We also proposed a particular time-stepping algorithm, which allowed us to completely decoupled local problems in time. The stability has been proven, with numerical experiments in accordance with the theory. Preliminary results concerning the convergence of the methods regarding both time and space variables have been conducted.

Perspectives

We envision the following possible extensions and perspectives for the different topics we addressed in this thesis.

Material models

Some of the simulation of direct bandgap materials in this work presented an absorption above the bandgap, which is unphysical. Indeed, this region is characterized by an exponential decrease of the imaginary part of the permittivity towards zero called exponential Urbach tail. This characteristic is not easily taken into account in the generalized dispersion model. Indeed, we showed some improvement in this region with very fine tuning of the simulated annealing algorithm, but there is still parasitic absorption. The modeling of such behavior requires the proposal of an other dispersion model, which can represent the Urbach tail. Exponential models seems to handle this features and could be interesting, but a preliminary study about its time-domain writing needs to be done.

Solution accuracy

The methodological contributions presented in this thesis were relatively simple, in the sense that they were based on simple criteria. The tuning of the different parameters generally relies on multiple experiments. Those strategies could be coupled with an a posteriori error estimate, which will allow an adaptation of the polynomial order as well as the refinement of the mesh in an automatic computation framework.

Photovoltaics and semiconductor physics

In this thesis, we treated the optical part of the photovoltaic. To simulate in a broader sense the absorption in solar cells, it would be interesting to consider, in addition to the optical model which characterizes the absorption of photons, an electrical model for the simulation of charge carriers transport. The electronic step can be represented threw a drift-diffusion electronic (DDE) model to calculate the current curve of the cell, and so its efficiency. This has been done in [ACM] by coupling a RCWA solver to a HDG solver.

Multiscale Hybrid-Mixed method

Several extensions are needed for the MHM method to compete and be fully compared to other numerical methods. At first, to deal with classical solar cell configurations, we should couple the method with the general dispersion model, and see how it is impacting both locals and global problems. In terms of performance, it would be interesting to work on different parallel versions of the method. Indeed, the divide-and-conquer framework of the method seemed to allow the use of several levels of parallelism: each local problems can be solved on its computer node, and

5

  Design optimization 5.1 Efficient Global Optimization (EGO) . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Design of Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Surrogate model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.3 Acquisition functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.4 EGO workflow and illustrative example . . . . . . . . . . . . . . . . . . . . 5.2 Surface texturing for enhanced color rendering . . . . . . . . . . . . . . . . . . . . . 5.2.1 Objectives of the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Cell composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3 Construction of geometrical models . . . . . . . . . . . . . . . . . . . . . . . 5.2.4 Influence of the geometry of the pyramids . . . . . . . . . . . . . . . . . . . 5.2.5 Optimization of the geometry of the pyramids . . . . . . . . . . . . . . . . . 5.3 Hot carrier pyramidal nanostructured solar cell . . . . . . . . . . . . . . . . . . . . 5.3.1 Material models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Comparison with FDTD results . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.3 Parallel EGO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.4 Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 High performance computing aspects 6.1 Coarse grain parallelism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Mesh partitioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.2 Parallel load balancing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.3 Iterative approach with advanced weighting . . . . . . . . . . . . . . . . . . 6.1.4 Strong scalability assessment . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.1 | General electromagnetic spectrum and zoom on the nanophotonics range of interest.

Figure 1 . 2 |

 12 Figure 1.2 | Solar PV generation and cumulative capacity by region, 2017-2023. Source: IEA website: https://www.iea.org/topics/renewables/solar/

Figure 1

 1 Figure 1.3 | Standard solar glass (left) versus light trapping (right). Source: Saint Gobain

Figure 2 .

 2 Figure 2.1 | Net power generating capacity added in 2017 by leading technology. Source: Frankfurt School-UNEP Centre and BNEF (2018) at SolarPower Europe 2018.

Figure 2 . 2 |

 22 Figure 2.2 | Components of a solar panel system Source: GreenMatch.co.uk

Figure 2 .

 2 Figure 2.3 | Diagram of the diffusion across a p-n junction, with the resultant uncovered space charges, the electric field and the drift currents. Source: Wikibooks.org

Figure 2 . 4 |

 24 Figure 2.4 | Energy vs. crystal momentum for a semiconductor with a direct bandgap (a) and an indirect bandgap (b).In the case where the energy of the photon is higher than the bandgap, an electron-hole pair can be generated by exciting an electron of the valence in the conduction band. Any extra energy from the photon above the bandgap is lost due to thermalization.

Figure 2 . 5 |

 25 Figure 2.5 | The absorption coefficient, for three different materials (GaAs, Si and aSi) at 300K depending on the wavelength of the light in vacuum.

Figure 2 . 7 |

 27 Figure 2.7 | Illustration of light trapping.

Figure 2 . 6 |

 26 Figure 2.6 | Efficiency of research solar cells classified by materials and types from 1975 to nowadays. Source: National Renewable Energy Laboratory (2018)

Figure 2 . 8 |

 28 Figure 2.8 | Real and imaginary parts of the microcrystalline silicon relative permittivity predicted by the generalized dispersive model compared to experimental data.

Figure 2 . 9 |

 29 Figure 2.9 | Linear mapping from the reference element T to the physical element Ti.

Figure 2 .

 2 Figure 2.13 | Ghost cells layer on the computational domain boundary. Boundary conditions are naturally handle via numerical fluxes and ghost cells.

Figure 2 .

 2 Figure 2.14 | General functioning of the PML.Outgoing wave propagates in the physical domain toward the PML, and crosses the interface. No reflection occurs, and the wave is progressively damped by the artificial medium and continues to propagate in the PML. At one point, the wave might encounters the edge of the computational domain, which is usually PEC or ABC. In both cases, the remaining of the wave will be totally (PEC) or partially (ABC) reflected toward the domain. In practice, the wave will travel a second time over the PML length, causing more damping. In general, when it enters again in the physical domain, the amplitude of the wave is attenuated by several orders of magnitude.

Figure 2 .

 2 Figure 2.15 | Wideband pulse representation in time-domain and frequency-domain. Parameters are E0 = 20 V.m -1 , ω0 = 3, t0 = 5. The shorter the time-domain pulse, the wider the frequency-domain spectrum.

  Locally adaptive DGTD method . . . . . . . . . . . . . . . . . . . . . . 54 3.1.1 Strategy based on the computed time-step per cell . . . . . . . . . . . . . 55 3.1.2 Strategy based on the number of points per wavelength . . . . . . . . . . 56 3.1.3 Performance assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 DGTD on hybrid meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.1 Hybrid configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.2 DGTD on hybrid meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.3 Numerical study of h-convergence . . . . . . . . . . . . . . . . . . . . . . 64 3.2.4 Numerical study of hybrid interfaces . . . . . . . . . . . . . . . . . . . . . 64 3.3 Realistic PV cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3.1 Numerical assessment for a double conical grating cell . . . . . . . . . . . 66 3.3.2 Numerical assessment for multilayer PV type nanostructures . . . . . . . 67 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 3

 3 Figure 3.1 | Nanostructured mesh of the double grating solar cell: the color corresponds to the cell size.

Figure 3 . 2 |

 32 Figure 3.2 | Comparison between the two locally adaptative strategies for P1 -P4 case. The repartition in the right figure shows the amount of cell per order, which is not directly related to hmin.

  Figure 3.3 | Reflection coefficient for DGTD-P4 (reference) and P1 -P4 strategies.

Figure 3 . 4 |

 34 Figure 3.4 | Non-conformity considered. The union of hybrid faces from the children cells needs to match perfectly the face of the parent neighbor.

Figure 3 . 5 |

 35 Figure 3.5 | Authorized hybrid interfaces.For the sake of precision, orders have to be higher on the parent element.

Figure 3

 3 Figure 3.6 | Element distribution for the h-convergence validation. To test convergence, the domain is arbitrarily cut in two halves, each part meshed with a different element type.

Figure 3 . 7 |

 37 Figure 3.7 | Hybrid configuration considered. Different hybrid interfaces to illustrate interface interferences.

4 Figure 3

 43 Figure 3.8 | Analysis of the error (L 2 in space, L ∞ in time) for different non conform interfaces (config 1 and config 3). Computation time Speedup P4 26 h 00 mn -P1 -P4, n ppw = 6 13 h 12 mn 1.95 P1 -P4 -Q3, n ppw = 6 4 h 26 mn 5.45 (3.00) P1 -P4 -Q4, n ppw = 6 9 h 12 mn 2.85 (1.40)Table 3.3 | Computation time for each simulation set-up on a Intel(R) Xeon(R) Gold 6154 CPU at 3.00GHz.

  Figure 3.9 | PV stack example with structured interface. The textured part represents a vary small percentage of the total domain.

Figure 3 .

 3 Figure 3.10 | Nanostructured mesh of the double grating solar cell: the color corresponds to the cell size.

Figure 3 .

 3 Figure 3.11 | Simulation of light trapping in a solar cell based on nanocone gratings. Simulations using hybrid mesh: contour lines of the module of the DFT of E for a wavelength λ = 450 nm.

Figure 3

 3 Figure 3.13 | Zoom on an interface region for the five considered meshes. Each color represents a specific material.

Figure 3 .

 3 Figure 3.14 | Convergence study for meshes M1 to M5 on a quantity of interest (A). Subfigure (f) compares the results obtained for each mesh.

Figure 4 .

 4 Figure 4.1 | A very coarse stratification (b) of a cone (a). Non-cartesian geometries in the z-direction are split in a set of z parallel layers.

  Yee cell (b) Staircasing effect on a complex interface

Figure 4 . 2 |

 42 Figure 4.2 | Yee cell and staircasing effect. (a) The H field components are on the center of the faces, while the E ones are on the center of the edges. (b) Cartesian discretization of a textured interface. The materials are spread in cells according to the percentage of volume physically above or below the textured surface.

Figure 4 . 3 |

 43 Figure 4.3 | The d.o.f. are associated to the edges of the element to naturally satisfy the continuity requirements for the approximated vector fields.

  Device enclosed between a TF/SF infinite plane and an infinite mirror

Figure 4 . 4 |

 44 Figure 4.4 | Periodic array of a textured layer.The TF/SF planes, Sin and Sout, are respectively above and below the layer (see the axes). The injection is done in the z direction and PML layers , while it is periodic in the x and y directions.

Figure 4

 4 Figure 4.5 | Standardized AM1.5G spectrum. This spectrum is a good representation of the illumination conditions in the geographical mid-latitudes on a tilted PV flat plate array on a clear day without clouds [GME02].

Figure 4 . 6 |

 46 Figure 4.6 | Convolution of the input topography data with a moving-average filter. Center element of the kernel is placed over the source data, which is replaced by a average of itself and the nearby values.

  .18) These filters have been assessed on a noisy set of points representative of PV random textures. We illustrate in figure 4.7 how Box blur strategy is way more invasive than the Gaussian smoothing. Of course, even if the blur box has succeeded in removing the unrealistic artefacts, it has also lost a large amount of geometric details, resulting in a very deep blur in figure 4.7(b). On the other hand, the Gaussian filter presents results much more consistent to the starting data. It allows us to remove the majority of defaults while keeping the texture details, as can be seen on the 3D topographic view in figure 4.8.

Figure 4 . 7 |

 47 Figure 4.7 | Comparison between the box blur and the Gaussian smoothing. The original image (a) presents some horizontal strips (artifacts caused by the microscope measures), (b) shows a box blur results, using a 5 × 5 kernel. The resulting image is blurry, which means a loss in texture details; Gaussian smoothing (c) is much closer to the original image with an attenuation of the strips.

Figure 4 .

 4 Figure 4.8 | Textured surface pre-and post-treatment Realistic texturation of solar cell layer.

  (a) Starting from a given unit cell (b) Merge with its reflection symmetry from x (c) Merge with its reflection symmetry from y

Figure 4 .

 4 Figure 4.9 | Symmetric mesh strategy.

Figure 4 .

 4 Figure 4.10 | Example of Tukey-window approach. As one can see in (b), the edges are all set-up to the mean value of a(x, y).

  (a) Original layer (b) Layer with PBC after applying Tukey window

Figure 4 .

 4 Figure 4.11 | Example of Tukey-window approach on a realistic textured layer with r = 0.3.

Figure 4 . 12 |

 412 Figure 4.12 | Composition of each structure. Material Structure 1 Structure 2 & 3 Structure 4 Structure 5 Air semi ∞ semi ∞ semi ∞ semi ∞ Si 3 N 4 0 60 60 60 AlGaAs 0 10 10 10 GaAs 100 100 100 100 AlGaAs 0 10 10 10 Al 2 O 3 0 10 10 10 TiO 2 (h) 100 100 100 210 Ag semi ∞ semi ∞ semi ∞ semi ∞

Figure 4 .

 4 Figure 4.13 | Mesh of the textured layer with a complexification: (a) corresponds to the symmetric pattern (w = 300 nm), (b) is breaking this symmetry (w = 450 nm), (c) is considering an evolution of the pattern in the z direction (w = 450 nm), and (d) is considering a pyramid (w = 450 nm).The red surface corresponds to PBCs, the gray is the interface between the Ag nanostructure and the substrate, and the green is the interface between the GaAs and the Ag.

Figure 4 .

 4 Figure 4.14 | Real and imaginary parts of the relative permittivity of Ag, GaAs and AlGaAs predicted by our dispersive model compared to experimental data.

Figure 4 .

 4 Figure 4.15 | Atot for Structure 1 using coarse mesh with polynomial orders ranging from 1 to 4.

Figure 4 .

 4 Figure 4.16 | Hierarchy of tetrahedral meshes for Structure 1.

Figure 4 .

 4 Figure 4.17 | Atot for Structure 1 using three gradually refined meshes with P4 interpolation.

Figure 4 .

 4 Figure 4.18 | A GaAs for Structure 1 using coarse mesh with polynomial orders ranging from 1 to 4.

Figure 4 .

 4 Figure 4.19 | A GaAS for Structure 1 using three gradually refined meshes with P4 interpolation.

Figure 4 .

 4 Figure 4.20 | Atot for the different structures with P3 and P4 interpolations on the coarse mesh.

Figure 4 .

 4 Figure 4.21 | A for all layers in the 5 different structures obtained using DGTD-P4 on a fine mesh.

Figure 4 .

 4 Figure 4.22 | Comparison of DGTD, FDTD and RCWA methods for the computation of the total absorption in the 5 structures.

Figure 4 .

 4 Figure 4.23 | Crystalline silicon structure with the cones extremity truncated.

Figure 4 .

 4 Figure 4.24 | Simulation using mesh M1: absorption spectra.

Figure 4 .

 4 Figure 4.25 | Simulation using mesh M2: absorption spectrum.

Figure 4 .

 4 Figure 4.26 | Comparaison of results using meshes M1 and M2: absorption spectra.

Figure 4 .

 4 Figure 4.27 | Simulation of light trapping in a solar cell based on nanocone gratings. Simulations using mesh M1: contour lines of the module of the DFT of E for a wavelength λ = 857 nm.

Figure 4 .

 4 Figure 4.28 | Simulation using mesh M2 and the DGTD-P2 method: contour lines of the module of the DFT of E for a wavelength λ = 857 nm.P Energy level max |A i+1 -A i | Physical time T Computation time #Iterations P 3 -P 4 5% / 3.21e-14 2 h 10 mn 2% 1.28e-02 4.52e-14 3 h 03 mn 1% 1.50e-03 6.02e-14 4 h 04 mn 0.5% 3.00e-03 1.00e-13 6 h 45 mn P 2 -P 3 5% / 3.82e-13 7 h 55 mn 2% 5.22e-01 (4.28e-01) 5.23e-13 10 h 50 mn 1% 3.36e-01 (5.82e-02) 6.52e-13 13 h 30 mn 0.5% 1.79e-01 (3.71e-02) 7.93e-13 16 h 25 mn 0.1% 3.94e-01 (4.92e-02) 1.20e-12 24 h 55 mn 0.05% 1.76e-01 (1.01e-02) 1.38e-12 28 h 32 mn 0.01% 1.90e-01 (8.25e-03) 1.93e-12 40 h 03 mn 0.005% 1.11e-01 (1.90e-03) 2.31e-12 47 h 56 mn 0.001% 8.69e-02 (1.46e-03) 3.10e-12 64 h 12 mn

Figure 4 .

 4 Figure 4.29 | Computation of the total absorption for different values of the threshold on the residual energy. We clearly see that in the lower part of the spectrum, the absorption converges faster than in the upper part.

Figure 4 .

 4 Figure 4.30 | Differences between two consecutive outputs of the absorption. One can see the two behaviors of the error in those graphs. In (a), three consecutive errors are plotted on the entire spectrum. In (b) and (c), we rescale the errors for the intervals [200 nm,650 nm] and [650 nm,1200 nm], respectively. There is three levels of magnitude between each.

Figure 4 .

 4 Figure 4.31 | Comparison of the obtained absorption with theoretical limits. A moving average with a window of 3 points is computed for easier readability. Theoretical limits are computed using optical parameters from literature.

Figure 4 .

 4 Figure 4.32 | Absorption obtained by two simulations and compared with the single run solution. One can see in (b) that the connection of the two simulations is almost perfect.

Figure 4 .

 4 Figure 4.33 | mesh of the GaAs cone with the injection region at the top and the PEC boundary at the bottom.

Figure 4 .

 4 Figure 4.34 | Zoom on cone extremity. Here the extremity of the cone was truncated. It is represented by a disk of radius r = 5nm.

Figure 4 .

 4 Figure 4.35 | Atot in the GaAs photonic crystal on coarse mesh with polynomial orders from 1 to 4.

Figure 4 .

 4 Figure 4.36 | Atot in the GaAs photonic crystal on fine mesh with polynomial orders from 1 to 4.

Figure 4 .

 4 Figure 4.37 | Original fitting of imaginary part of the GaAs electric permittivity.

Figure 4 .

 4 Figure 4.38 | Optimized fitting of imaginary part of the GaAs electric permittivity. We observe a clear improvement near the band gap region.

Figure 4 .

 4 Figure 4.39 | Fitting of imaginary part of the GaAs electric permittivity in semi-log scale. (Left) is the original fit while (Right) is the optimal fit.

Figure 4 .

 4 Figure 4.40 | Simulation using DGTD-P2 method. (Left) is based on the original fit while (Right) is based on the optimal fit. The obtained non physical absorption is divided by almost 3.

Figure 4 .

 4 Figure 4.41 | Simulations using DGTD-P2 method and second option for the fitting: contour lines of the module of the DFT of E (top) and real part of Ex (bottom) for a wavelength λ = 800 nm.

Figure 4 .

 4 Figure 4.42 | Composition of the thin-film solar cell device.

  (a) Modified mesh using Tukey window approach (b) Modified mesh using symmetric mesh

Figure 4 .

 4 Figure 4.43 | Submodel meshes obtained with the two different strategies for imposing periodicity.

Figure 4 .

 4 Figure 4.44 | Full cell using the Tukey window approach with r = 0.3.

  µc-Si:H, imaginary part

Figure 4 .

 4 Figure 4.45 | Real and imaginary parts of the relative permittivity of front TCO, a-Si:H, µc-Si:H and back TCO predicted by our dispersive model compared to experimental data.

Figure 4 .

 4 Figure 4.46 | EQE for 1 × 1 µm 2 flat submodel (mesh M1) with polynomial orders from 1 to 4.

Figure 4 .

 4 Figure 4.47 | EQE for 1 × 1 µm 2 submodel (mesh M2) with polynomial orders from 1 to 4.

Figure 4 .

 4 Figure 4.48 | EQE for 10 × 10 µm 2 model (mesh M3) with polynomial orders from 1 to 4.

Figure 4 .

 4 Figure 4.49 | Comparison of EQE from simulations with the FDTD and DGTD solvers, for microcrystalline silicon (µc-Si:H) and amorphous silicon (a-Si:H) using the full model (mesh M2).

Contents 5 . 1

 51 Efficient Global Optimization (EGO) . . . . . . . . . . . . . . . . . . . 117 5.1.1 Design of Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.1.2 Surrogate model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.1.3 Acquisition functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.1.4 EGO workflow and illustrative example . . . . . . . . . . . . . . . . . . . 118 5.2 Surface texturing for enhanced color rendering . . . . . . . . . . . . . 122 5.2.1 Objectives of the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.2.2 Cell composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.2.3 Construction of geometrical models . . . . . . . . . . . . . . . . . . . . . . 124 5.2.4 Influence of the geometry of the pyramids . . . . . . . . . . . . . . . . . . 126 5.2.5 Optimization of the geometry of the pyramids . . . . . . . . . . . . . . . . 128 5.3 Hot carrier pyramidal nanostructured solar cell . . . . . . . . . . . . . 131 5.3.1 Material models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.3.2 Comparison with FDTD results . . . . . . . . . . . . . . . . . . . . . . . . 133 5.3.3 Parallel EGO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.3.4 Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

Figure 5 . 1 |

 51 Figure 5.1 | Illustration of LHS sampling for 2 discrete parameters with n = 20.

Stop•Figure 5 . 2 |

 52 Figure 5.2 | EGO optimization flow diagram. The orange part represents the iterative part of the EGO.

Figure 5

 5 Figure 5.3 | Simple 1D example of the EGO. Superposition of the analytical objective function with the surrogate for several iterations of the EGO algorithm (Left). Plot of the expected improvement in orange and of the future sampling point (Right).

  Distance between consecutive points.

Figure 5 . 4 |

 54 Figure 5.4 | Simple 1D example of the EGO.In the left we can see the distance between two consecutive x + . The iterations 1 to 5 correspond to the exploration part. In the right we plot the current maximum.

Figure 5

 5 Figure 5.5 | Illustration of the problem of poor absorption of red tints (Left). Examples of texturations of a PV glass structure (Right) -Gray: back refelector (metal) -Brown: TCO (Transparent Conducting Oxide) -Green: glass -Dark blue: anti-reflective gel coating -Light blue: ambient medium.

Figure 5

 5 Figure 5.6 | Structure of a classical solar cell stack: top and transverse views.

Figure 5 . 7 |

 57 Figure 5.7 | Several nanostructurations of a ZnO layer surface (top view).

  Figure 5.8 | Composition of the PV cell stack.

  .10(a). It is also possible to reduce Material 1 | Layers thicknesses. Glass thickness is 500 µm.

Figure 5 . 9 |

 59 Figure 5.9 | Geometrical model of the PV stack with a texturing of the layers based on a pyramidal grating. This corresponds to a 5 × 4 pyramidal array. For the computation, we consider one pyramid with PBCs.

  Figure 5.10 | Unit pattern with an isolation of the very thin Si-based layers.

Figure 5 .

 5 Figure 5.11 | Images of the nanotexturing of a AZO layer surface (top view) used in the definition of the selected texture.

Figure 5 . 12 |

 512 Figure 5.12 | Absorption coefficient convergence for the starting-point texture with h=47.8 nm and r=165.1 nm.

Figure 5 .

 5 Figure 5.13 | Absorption coefficient for various textures obtained with DGTD-P3 method. Computation times for each configuration are respectively 1181 sec, 1360 sec and 1356 sec.

  Figure 5.14 | Absorption coefficient for various textures with DGTD-P3 method. Computation times for each configuration are respectively 1181 sec, 1693 sec, 2017 sec.

  V3: h=175.8 nm, r=192.3 nm.

Figure 5 .

 5 Figure 5.16 | Absorption coefficient for the starting-point texture with h=47.8 nm and r=165.1 nm, and the 3 optimized textures for DGTD-P3 method.

Figure 5 .

 5 Figure5.17 | Simulations using DGTD-P3 method on the optimal configuration: contour lines of the real part of Ex for a wavelength λ = 850 nm (top) and λ = 400 nm (bottom).

2 Figure 5

 25 Figure 5.18 | Composition of the PV cell. Gray area corresponds to artificial volumes (PML and SF region).Orange area corresponds to both the absorber (InGaAs layer) and the nanotexture (SiO2 pyramid).

Figure 5 . 19 |

 519 Figure 5.19 | Illustration of the three design parameters. The gap is the smallest distance between two successive pyramids.

Figure 5

 5 Figure 5.20 | Real and imaginary parts of the relative permittivity of Ag, InGaAs, MgF2 and ZnS predicted by our dispersive model compared to experimental data.

  Figure 5.21 | Real and imaginary parts of the relative permittivity predicted by our dispersive model compared to experimental data in semi-log scale.

Figure 5

 5 Figure 5.22 | Convergence study on configurations D1 and D2 compared with FDTD.

Figure 5

 5 Figure 5.23 | Superposition of our results with theoritical limits.

Figure 5 .

 5 Figure 5.24 | q-EGO optimization flow diagram. The orange part represent the iterative iterative portion of the EGO method.

  Figure 5.25 | Prediction of the Jsc from the surrogate model after DoE step and during several iterations of the parallel EGO algorithm. Squares represents the DoE samples and the circles are the point selectionned by the algorithm.

Figure 5

 5 Figure 5.26 | Comparison of the initial design and the optimal one.

Figure 5 .Figure 5

 55 Figure 5.27 | Simulations using DGTD-P4 method on the optimal configuration: contour lines of the real part of Ex (bottom) for a wavelength λ = 1290 nm.

Contents 6 . 1

 61 Coarse grain parallelism . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.1.1 Mesh partitioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.1.2 Parallel load balancing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.1.3 Iterative approach with advanced weighting . . . . . . . . . . . . . . . . . 148 6.1.4 Strong scalability assessment . . . . . . . . . . . . . . . . . . . . . . . . . 151 6.2 Hybrid parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 6.2.1 Performance results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Figure 6 . 1 |

 61 Figure 6.1 | Mesh partitioning workflow. Starting from a mesh (a), its corresponding dual graph is computed.Then some partitioning algorithms are used to obtain partitions of this dual graph (c). From here, the submeshes are recovered.

Figure 6 . 4 |

 64 Figure 6.4 | Tetrahedral (left) and hybrid cubic-tetrahedral (right) meshes of a nanostructured solar cell.

Figure 6

 6 Figure6.5 | Weak scalability comparison between a DGTD simulation using the hybrid mesh M1 and the tetrahedral mesh M2. The sequential simulation is 60% faster with the hybrid mesh, however, the parallel speedup is much worse to the point of being less interesting when the number of subdomains exceeds 4.

Figure 6 . 7 |

 67 Figure 6.7 | Workflow of the iterative partitioning method.

  Figure 6.8 | Illustration of the iterative partition steps. The gradient of blue represent the weight associate to the underlying element.

Figure 6 .

 6 Figure 6.11 | DIOGENeS-DGTD traces: focus on the main time loop. Isolation of one iteration of DIOGENeS-DGTD with LRSK-14-5 time scheme on 8 cores. The highlighted iteration corresponds to the update of the reflection observable.

Figure 7 . 1 |

 71 Figure 7.1 | Space discretization. First level space discretization. The domain Ω is split into coarse macroelements K.

Figure 7 . 2 |

 72 Figure 7.2 | Time discretization. First level time discretization. The time interval [0, T ] is split in several timeslab In. The blue H is the maximum size of K1, the orange H is the maximum size of K2.

Figure 7 .

 7 Figure 7.3 | Multiscale basis function. Example of η n,E i basis function on a given macroelement K.

Figure 7 . 5 |

 75 Figure 7.5 | Time discretization. Two level time discretization. The timeslab ∆t is split in micro time-steps ∆τ

4 :if n = 1 then 6 : 7 :

 467 for m = 1 → M do Loop on local time-steps 5:Determine η m,H from η m-1 2 ,E and η m-1,H Determine η m+ 1 2 ,E from η m-1 2 ,E and η m,H

p e = 4 Figure 7 . 6 |

 476 Figure 7.6 | Convergence study of MHM-DGTD method for all the order couples (p e , p K ) with p e from 1 (a) to 4 (d) and p K from 1 (orange) to 4 (light orange).

Figure 7 . 7 |

 77 Figure 7.7 | Energy preservation of MHM-DGTD method with LTS and two different integration in time for ∆τ constant and increasing ∆t.

Figure 7 . 8 |

 78 Figure 7.8 | Convergence study of MHM-DGTD method with LTS with (p e , p K ) = (1, 3), for fixed "small enough" (h, ∆τ, H) with H = 3.12e-2, ∆τ = Ch 10 and ∆t → 0 (right to left).

Figure 7 . 9 |

 79 Figure 7.9 | Convergence study of MHM-DGTD method with LTS with (p e , p K ) = (1, 3) with (∆t, ∆τ, H, h) → 0 for spatial refinement from 1 to 8 (h = H α for α ∈ [1, 8]).

Figure 7 .

 7 Figure 7.10 | Convergence study of MHM-DGTD method with LTS for fixed small enough local problems K(h, ∆τ ) with p e from 1 (a) to 4 (d) and p K from 1 (orange) to 4 (light orange).

7. 4 . 4

 44 Comparison with DGTD 7.4.4.1 Photonic structure

Figure 7 .

 7 Figure 7.11 | Geometry of the nanowaveguide considered. The surrounding is made of air (orange), the chip is composed of silica (blue) with holes of silicium (dark blue).

Figure 7 . 12 |

 712 Figure 7.12 | Points on which the errors is computed. We selected 9 points uniformly distributed in order to compute the error between a reference solution and the bench runs.

Figure 7 .

 7 Figure 7.13 | Meshes used for the methods comparison. All meshes are aligned with the silicium device but M1 is not conform with the spheres.

Figure 7 .

 7 Figure 7.15 | Representation of the computational times evolution when increasing ∆t with the MTS time scheme.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 8.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 8.2.1 Material models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 8.2.2 Solution accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 8.2.3 Photovoltaics and semiconductor physics . . . . . . . . . . . . . . . . . . 187 8.2.4 Multiscale Hybrid-Mixed method . . . . . . . . . . . . . . . . . . . . . . . 187 Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197

  H and A α,E are equal, e.g. A α,H = A α,E ≡ A α ; 2. A H and A E are multiples of each other, e.g. A H = εr µr A E ≡ Y 2 εr µr A U .

  is the estimate obtained by the forward Euler method. The slope of the solution between tn and tn + ∆t is approached by φ1. The RK2 estimate is obtained by averaging φ1 with a second approximation of the slope obtained from φ n+1 * in tn + ∆t. φ represents the average of the different slopes.

	φ n+1
	φ(tn+1)	φ n+1				φ ( t )
	φ n+1 *	φ n φ 1	tn	φ φ2	t n+ 1 2 φ2	2 φ ( t ) φ φ 2	tn+1
	φ n			φ1
	tn Figure 2.12 | Steps of the LSRK2 algorithm. Here, φ2 represents ∆tf tn + ∆t tn+1 2 , φ 1 with 2 stages.
	Figure 2.10 | Steps of the RK2 algorithm. φ n+1 *
						φ 4
	φ n+1
					φ	φ ( t )
				φ 3
	φ n		φ1	φ 2
	tn			t n+ 1 2		tn+1
	Figure 2.11 | Steps of the RK4 algorithm. φ represents the average of the different slopes.

4 :

 4 Sort cells by ascending ∆t i 5: p(1) ← p min 6: ∆t loc ← ∆t p min 1 7: p loc ← p min 8: for i ← 2, n t do Go over cells by ascending order of ∆t i

	9:	if p loc + 1 > p max then	Check that we do not exceed p max
	10: 11: 12:	p(i) ← p max else Compute ∆t p loc +1

i 13: if ∆t p loc +1 i > ∆t loc then Check if it is worth changing order 14: ∆t loc ← ∆t p loc +1 i Update the limiting time-step 15: p(i) ← p loc + 1 16: p loc ← p loc + 1 Update the current order 17: else 18: p(i) ← p loc 19: end if 20:

  28 × 10 -2 1.98 3.27 × 10 -3 1.97 1.47 × 10 -3 2.03 P1 -Q2 6.90 × 10 -2 -1.88 × 10 -2 1.87 4.82 × 10 -3 1.96 2.16 × 10 -3 2.05 P1 -Q3 6.99 × 10 -2 -1.89 × 10 -2 1.89 4.83 × 10 -3 1.97 2.16 × 10 -3 2.06 P1 -Q4 6.69 × 10 -2 -1.89 × 10 -2 1.89 4.83 × 10 -3 1.97 2.16 × 10 -3 2.06 -4.02 × 10 -5 4.04 2.54 × 10 -6 3.99 5.13 × 10 -7 4.07 P3 -Q4 6.34 × 10 -4 -4.07 × 10 -5 3.96 2.56 × 10 -6 3.99 5.18 × 10 -7 4.08 Table 3.2 | Error levels and convergence rates of the cubic cavity case for mixed orders of approximation on hybrid meshes of increasing refinement. All simulations were run with upwind fluxes and LSRK4 time integration.

	P2 -Q2 P2 -Q3 P2 -Q4 P3 -Q3	7.26 × 10 -3 -9.82 × 10 -4 2.89 1.26 × 10 -4 2.97 3.74 × 10 -5 3.09 7.07 × 10 -3 -9.59 × 10 -4 2.88 1.23 × 10 -4 2.96 3.68 × 10 -5 3.08 7.04 × 10 -3 -9.67 × 10 -4 2.86 1.26 × 10 -4 2.94 3.79 × 10 -5 3.05 6.62 × 10 -4

table 3 .

 3 4).

		Number of cells hmin (nm) hmax (nm)
	M1 M2 M3 M4 M5	24000 hexahedra 320000 hexahedra 51816 tetrahedra 414528 tetrahedra 5094 hexahedra + 21380 tetrahedra	40 20 11 6 6	70 20 25 19 100

Table 3 .

 3 4 | Characteristics of different configuration meshes.

  4.6 Thin-film solar cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.6.1 Geometrical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.6.2 Material models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.6.3 Numerical convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.6.4 Comparison with FDTD simulations . . . . . . . . . . . . . . . . . . . . . 112 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Table 4 .

 4 1 | Thicknesses of the different materials for Structure 1 to 5 in nanometers.

	Mesh Structure 1 (very coarse) Structure 1 (coarse) Structure 1 (fine) Structure 2 Structure 3 Structure 4 Structure 5	Tetrahedron h min (nm) h max (nm) hmax h min 314 100.0 327.7 3.28 3638 31.5 164.5 5.23 75135 13.9 57.2 4.11 95052 10.0 54.4 5.43 95062 10.0 54.8 5.48 108102 10.0 57.1 5.71 106656 10.0 79.0 7.91

Table 4 .

 4 2 | Characteristics of tetrahedral meshes of Structures 1 to 5.

Table 4 .

 4

	hmax
	h min

4 | Characteristics of tetrahedral meshes of the conical nanowire .

Table 4

 4 .47, while the EQE profiles in figure4.48 are smoothed. From those statements, we can say that flat interfaces cause those resonances. It is also visible on the submodel results because of the Tukey window, which is flattening a large portion of the surface.

	Mesh M1: 1 × 1 µm 2 M2: 1 × 1 µm 2 M3: 10 × 10 µm 2	hmax h min 15.9 48.7 917.8 137.0 # Tetrahedra h min (nm) h max (nm) 29461 29.2 465.7 41387 9.9 482.6 1151793 6.7

.5 | Characteristics of the tetrahedral meshes for the 1 × 1 µm 2 and 10 × 10 µm 2 models.

  allows to perfectly sample the spectrum of the A layer operator, by evaluating equation (4.11) at certain ∆t obs time-steps. By doing so, we can reduce the CPU time up to 400% for the full model as one can see in table 4.6. These simulations have been performed on an in-house cluster system with 128 cores (16 Intel E5-2670 @ 2.6 GHz nodes each with 8 cores). Table 4.6 | CPU times with and without applying Shannon theorem (full model, mesh M2).

	Mesh	Order	Improved time	Original time Gain (%)
	1 × 1 µm 2	P 1 P 2 P 3 P 4	12 mn 38 s 19 mn 01 s 46 mn 58 s 3 h 25 mn 34 s	13 mn 58 s 23 mn 39 s 1 h 01 mn 08 s 3 h 45 mn 06 s	10.6 24.4 30.2 9.5
	10 × 10 µm 2 P 1 P 2 P 3	3 h 02 mn 16 s 9 h 29 mn 52 s 32 h 43 mn 57 s 126 h 12 mn 00 s 15 h 54 mn 14 s 49 h 46 mn 51 s	423.5 424.1 285.5

Table 5 .

 5 2 | Texturation of the AZO layer: geometrical characteristics of pyramids.

  a linear speed-up is observed (see table 6.2). It is interesting to see how the computation of this figure of merit makes the computation way more expensive. Table 6.1 | Strong scalability analysis of the DGTD-P k solver on the Occigen system. Tetrahedral mesh with 305,265 vertices and 1,689,764 elements. Timings for 1000 time iterations including the evaluation of the volume absorption in the Silicon layer. Execution mode: 1 MPI process per core. Table 6.2 | Strong scalability analysis of the DGTD-P k solver on the Occigen system. Mesh M1 (full model) with 305,265 vertices and 1,689,764 elements. Timings for 1000 time iterations excluding the evaluation of the volume absorption in the Silicon layer. Execution mode: 1 MPI process per core. Table 6.3 | Strong scalability analysis of the DGTD-P k solver with an optimized partition. Tetrahedral mesh with 305,265 vertices and 1,689,764 elements. Timings for 1000 time iterations including the evaluation of the volume absorption in the Silicon layer. Execution mode: 1 MPI process per core.

	Solver	# cores	Elapsed	Speedup
	DGTD-P1	96 192 384	2681 sec 1.00 (1.0) 1365 sec 1.95 (2.0) 768 sec 3.50 (4.0)
	DGTD-P2	96 192 384	4364 sec 1.00 (1.0) 2254 sec 1.95 (2.0) 1332 sec 3.30 (4.0)
	DGTD-P3	192 384	3678 sec 1.00 (1.0) 2232 sec 1.65 (2.0)
	Solver	# cores	Elapsed	Speedup
	DGTD-P1	96 192 384	584 sec 1.00 (1.0) 292 sec 2.00 (2.0) 146 sec 4.00 (4.0)
	DGTD-P2	96 192 384	974 sec 1.00 (1.0) 490 sec 2.00 (2.0) 246 sec 3.95 (4.0)
	DGTD-P3	192 384	808 sec 1.00 (1.0) 418 sec 1.95 (2.0)
	Solver	# cores	Elapsed	Speedup
	DGTD-P1	96 192 384	2217 sec 1.00 (1.0) 1114 sec 2.00 (2.0) 559 sec 3.95 (4.0)
	DGTD-P2	96 192 384	3836 sec 1.00 (1.0) 1880 sec 1.95 (2.0) 937 sec 3.90 (4.0)

  Table 6.5 | Performance results for the solar cell device test problem.

	Method	# time-steps # OpenMP tasks Wall time (sec) Speedup
	DGTD-P 1 -DGTD-P 2 ---DGTD-P 3 --DGTD-P 4 ---	80 -110 ---138 --179 ---	1 4 1 4 8 16 4 8 16 4 8 16 32	412.0 107.0 936.0 235.5 122.1 69.5 483.8 253.2 143.0 1124.5 568.1 323.5 221.0	1.00 3.85 1.00 3.98 7.67 13.47 1.00 1.91 3.38 1.00 1.98 3.48 5.10

  7.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183

  Let us define the second local level energy, ∀n ∈ {1, ..., N }, ∀m ∈ {1, ..., M K }

		, µ H ∂T H	= 0	(7.30)
			ˆIn
	which corresponds to a trapezoidal quadrature rule on the integral	E H,h , µ H ∂T H	= 0. This
	choice comes directly from the theorem below.
	Theorem 7.3.1. E	n,m+1/2 H,h

  We obtain the following result. We start from the fully discrete local energy, and by using the second equation of (7.36) to replace H n,m+1

	Theorem 7.3.2. For all K ∈ T H , the local energy E H,h n,m+1/2 H n,m H,h under the following CFL condition per macro element K (K) is a positive form of E n,m+ 1 2 H,h	and
					∆τ K <	h K 2C K	(7.40)
	H,h	, we get:
		E H,h n,m+1/2	(K) = =	ε 2 ε 2	E E	n,m+ 1 2 H,h n,m+ 1 2 H,h	2 0,T K h 2 0,T K h	+ +	1 2 µ 2	µH n,m H,h , H n,m+1 H,h H n,m H,h 2 0,T K h	h T K
				-∆τ K E	n,m+ 1 2 H,h	, ∇ × H n,m H,h T K h
				+∆τ K E n,m+ 1 2 H,h	, H n,m H,h	F K h

with C K a given constant independent of (∆τ K , h K );

Proof.

  Then, for all macro element K, the discrete electromagnetic energy E

								m+ 1 2 H,h	2 0,τ + µ H n,m H,h	2 0,τ
			= =	1 2 1 2	--	∆τ K h K ∆τ K h K	( C 2 C1 + C2 + C3 4ζ K h ζ K τ ∈T K ε E n,m+ 1 ) τ ∈T K h 2 H,h 2 0,τ + µ H n,m ε E n,m+ 1 2 H,h H,h	2 0,τ + µ H n,m H,h 2 0,τ	2 0,τ
	form of E	n,m+ 1 2 H,h	and H n,m H,h under the condition:	n,m+1/2 H,h	(K) is a positive
								∆τ K <	h K 2C	ζ K	(7.41)
	Therefore, the global energy E n H,h =	K∈T H	M K m=1	E	n,m+1/2 H,h

(K) is a positive form of all variables under a local CFL condition for each K ∈ T H .

Table 7 .

 7 1 | Comparison between MHM and DGTD on several meshes. Mesh Number of d.o.f. max |E -E ref |

	MHM DGTD-P4	M1 M1 M2 M3 M4	2450 7230 21690 107730 184260	4.82 × 10 -2 2.27 8.82 × 10 -1 2.79 × 10 -1 5.15 × 10 -2

The International Energy Agency

Institut Photovoltaïque d'Île-de-France

https://refractiveindex.info/

http://gmsh.info/

http://www.meshgems.com/

http://gmsh.info/

http://famosa.gforge.inria.fr/

https://cran.r-project.org/web/packages/DiceOptim/index.html

http://gmsh.info/

https://www.lumerical.com/products/fdtd/

https://diogenes.inria.fr/

https://tools.bsc.es/extrae

https://tools.bsc.es/paraver

Main supercomputer in the BSC

present in an iteration of the LSRK 14-5 scheme, presented in [START_REF] Niegemann | Efficient low-storage runge-kutta schemes with optimized stability regions[END_REF]. The iteration considered is an iteration for which the observed quantity is updated. As explained in chapter 4, this quantity is computed on the TF/SF interface, i.e. in the 6th subdomain. Since R is a frequency-dependent quantity, it does not need to be updated at each time-step according to Shannon-Nyquist sampling theorem, which states that a broadband time signal can be sampled and entirely reconstructed from its samples if the wave is sampled at least twice as fast as its highest frequency component [START_REF] Shannon | Communication in the presence of noise[END_REF]. We purposely removed the communications between cores while keeping the gap between each step, which is representative of the idle time. We clearly see the synchronization at each step. We observe an imbalance in the computational time for the 6th CPU, at the last step. This is relevant to the placement of the TF/SF interface, knowing that the observable is updated on this surface in the mesh. Furthermore, we see that the iteration responsible for the observable updates (blue) is even more imbalanced.

Those load balancing issues are even more problematic when dealing with hybrid meshes. We remind some results from chapter 3, where we have shown a gain of almost 40% for the -10

.9 | Correction of the load imbalance using the iterative partitioning approach. Each bar corresponds to the relative imbalance of a single CPU to the average value computed over all processors. As can be seen, the balance of the CPU loads is restore over iterations, with at the end a maximum imbalance of 3%.

time iterations. We consider the case of 8 partitions. For a simpler comparison of the different partitions, we are using the relative deviation to the mean CPU time ∆ (in %) on each process. The effect of weighting is assessed by analyzing the evolution of the relative deviations from the first iteration to the tenth. The results are shown in figure 6.9. As can be seen, the load balancing is improving (from 15% to 2%).

Then, we have used this strategy on the case shown in section 4.5.1, which exhibits notable load imbalance caused by the TF/SF interface included in one subdomain. Using this strategy over 20 time iterations allows us to restore an acceptable level of load imbalance. As one can see in figure 6.10, the new partition is splitting the TF/SF interface. The function traces in figure 6.11 highlight the better overall load balancing, even if the two aforementioned partitions are still the most expensive ones. An improvement of 10% is observed for this problem. boundary condition on ∂Ω given by:

Numerical study of the MTS scheme

Stability We illustrate the preservation of energy (with J = 0) as proved in Theorem 7.3.1. We show in figure 7.7 that energy is numerically preserved when a trapezoidal quadrature rule is used to evaluate the global problem. We also tested the choice of rectangular quadrature to approximate the global problem and even did not observe stability (see orange curve in figure 7.3.1).

Numerical convergence Since the theoretical analysis of convergence has not yet been done, the majority of this part consists of assumptions, factual observations and postulates. From classical analysis, one could naturally expect an error estimate as

The expectation concerning the H p e +1 comes from [START_REF] Harder | A hybrid-mixed method for elasticity[END_REF]. Our first intuition, especially following the degenerated tests and the use of the LF2 time scheme, is that min(k, e) = 2. The following tests aim at showing the behavior of the error when using the MTS time scheme, i.e. ∆t = ∆τ K for at least one K.

• Error analysis with respect to ∆t (figure 7.8) First, we choose (H, h K , ∆τ K ) for all K ∈ T H , small enough so as to decouple the error analysis and provide an error analysis with respect to ∆t. To do this, we start by choosing the pair of orders (p e , p K ) = (1, 3) such that the local problems are well posed following [START_REF] Harder | A hybrid-mixed method for elasticity[END_REF] (p k = p e + 2) without refinement at the second level, i.e. H = h K , ∀K ∈ T H . We also set ∆τ as one tenth of the optimal micro time-step which is defined as min K∈T H CF L h K . We start from the degenerated case i.e. ∆t = ∆τ K , ∀K ∈ T H , and then we gradually increase the ∆t. It was done for two different levels of spatial refinement, i.e. H/h K = 1 for the orange plot, and H/h K = 2 for the blue one. A second-order slope is observed for the bigger timesteps, followed by a saturation of a magnitude close to H 2 . This first observation confirms our intuitions that we also have a second order in ∆t i.e. e = 2.