A Robust GPU-Accelerated Kernel Compensation Solver with Novel Discretization for Photonic Crystals in Anisotropic Media

1. Énoncé du problème

L'article aborde le défi computationnel consistant à résoudre le problème aux valeurs propres de Maxwell pour des cristaux photoniques (CP) 3D composés de milieux anisotropes.

• Équations gouvernantes : Le problème consiste à trouver les fréquences propres ($\omega$) et les champs magnétiques ($\mathbf{H}$) satisfaisant $\nabla \times (\varepsilon^{-1} \nabla \times \mathbf{H}) = \omega^2 \mathbf{H}$ avec une condition de divergence nulle ($\nabla \cdot \mathbf{H} = 0$) et des conditions aux limites quasi-périodiques définies par le théorème de Bloch.
• Défis clés :
  1. Espace nul : L'opérateur rotationnel-rotationnel possède un espace nul massif (champs gradients), ce qui complique les solveurs de valeurs propres.
  2. Anisotropie : Lorsque le tenseur de permittivité ($\varepsilon$) contient des entrées hors-diagonale non nulles, les schémas aux différences finies standards font face à un désalignement des degrés de liberté (DoF) sur des grilles décalées (schéma de Yee), rendant difficile la construction d'un système défini positif hermitien (HPD).
  3. Évolutivité : La résolution de problèmes 3D à grande échelle nécessite une gestion efficace des matrices clairsemées et une parallélisation massive.

2. Méthodologie

Les auteurs proposent un cadre complet combinant le schéma aux différences finies de Yee, la compensation de noyau et l'accélération GPU sans matrice.

A. Discrétisation et transformation de coordonnées

• Schéma de Yee : Le champ magnétique est discrétisé sur les DoF des faces, et le champ électrique sur les DoF des arêtes.
• Transformation de coordonnées : Pour gérer des vecteurs de réseau arbitraires (SC, FCC, BCC), une transformation de coordonnées mappe le domaine physique sur un cube unité $(0,1)^3$, alignant la périodicité avec la base standard. Cela préserve la structure circulante des opérateurs différentiels.
• Opérateurs sans matrice : Au lieu de stocker des matrices clairsemées, les opérateurs différentiels sont implémentés en utilisant la Transformée de Fourier Discrète 3D (DFT). Cela réduit le coût de la multiplication matrice-vecteur à une seule DFT globale 3D par itération.

B. Technique de compensation de noyau

Pour éliminer l'espace nul de l'opérateur rotationnel tout en imposant la condition de divergence nulle, les auteurs emploient une formulation pénalisée :
$$(\mathcal{A}\mathcal{M}_\varepsilon\mathcal{A}^\dagger + \gamma \mathcal{B}^\dagger\mathcal{B})\mathbf{H}_h = \omega_h^2 \mathbf{H}_h$$

• $\mathcal{A}$ : Opérateur rotationnel discret.
• $\mathcal{B}$ : Opérateur divergence discret.
• $\gamma$ : Un coefficient de pénalité.
• Contribution théorique : L'article prouve rigoureusement que pour un $\gamma$ suffisamment grand, les plus petites valeurs propres positives de ce système pénalisé coïncident avec celles du problème original. Il fournit une règle explicite pour choisir $\gamma$ en fonction du vecteur d'onde $\mathbf{k}$ et de la taille de la grille, assurant que le système reste bien conditionné (contrairement aux méthodes FEM qui nécessitent souvent $\gamma \sim O(h^{-3})$).

C. Nouvelle discrétisation pour les milieux anisotropes

L'innovation centrale réside dans la gestion des entrées hors-diagonale du tenseur de permittivité $\varepsilon_1$ sur une grille décalée :

1. Le problème : Les termes hors-diagonale (par exemple $\varepsilon_{12}E_2$) couplent des composantes de champ définies sur des emplacements de grille différents (arêtes vs faces), entraînant un désalignement des DoF.
2. Solutions proposées :
   • Méthode triviale : Ignore le décalage spatial (approxime les DoF comme alignés).
   • Interpolation Cross-DoF (Nouvelle) : Utilise des opérateurs d'interpolation symétriques pour transférer les valeurs entre les grilles décalées. Crucial, les auteurs construisent un opérateur symétrisé pour garantir que la matrice discrète résultante $\mathcal{M}_\varepsilon$ reste hermitienne.
3. Preuve HPD : Les auteurs prouvent que sous des contraintes physiques (valeurs propres de $\varepsilon_1 \in (0, 1]$ et dominance diagonale spécifique ou conditions d'entrées nulles), la matrice de permittivité discrète résultante est définie positive hermitienne (HPD). Ceci est essentiel pour la stabilité de la méthode de compensation de noyau.

D. Extension d'ordre élevé et préconditionnement

• Maillages d'ordre élevé : Le cadre prend en charge des maillages aux différences finies d'ordre pair arbitraire (jusqu'à l'ordre 10 testé) sans briser la structure sans matrice.
• Préconditionnement : Un préconditionneur diagonal est dérivé dans le domaine de Fourier, permettant une inversion triviale et assurant que le nombre de conditionnement du système est borné par le nombre de conditionnement de la matrice de permittivité.

3. Contributions clés

1. Théorie rigoureuse de la compensation de noyau : Fournit une fondation théorique complète pour la méthode de compensation de noyau, incluant une caractérisation précise de l'espace nul et une formule explicite et indépendante de la grille pour le coefficient de pénalité $\gamma$.
2. Discrétisation anisotrope novatrice : Introduit une méthode mathématiquement prouvée pour discrétiser les tenseurs de permittivité hors-diagonale sur les grilles de Yee tout en préservant la propriété HPD, résolvant un problème de longue date dans les simulations de cristaux photoniques aux différences finies.
3. Accélération GPU sans matrice : Implémente l'intégralité du solveur (opérateurs, préconditionnement et produits matrice-vecteur) en utilisant la DFT 3D sur GPU, éliminant le besoin de stockage de matrices clairsemées et permettant des accélérations massives.
4. Précision d'ordre élevé : Étend le cadre sans matrice à des discrétisations d'ordre élevé, améliorant la précision pour les problèmes lisses sans augmenter significativement la surcharge mémoire.

4. Résultats numériques

La méthode a été validée sur diverses structures de réseau (SC, FCC, BCC) avec des interfaces courbes et gyroïdes, à la fois isotropes et anisotropes.

• Précision :
  • La méthode d'interpolation "Cross-DoF" produit une précision légèrement supérieure à la méthode "Triviale".
  • Pour des distributions de permittivité lisses, le maillage d'ordre 10 a atteint une convergence d'environ l'ordre 10, confirmant les attentes théoriques.
  • Les rapports de bandes interdites et les points de Weyl ont été reproduits avec précision, correspondant à la littérature existante.
• Performance :
  • Accélération : L'accélération GPU a atteint une accélération moyenne de >40x par rapport aux implémentations CPU (en utilisant NumPy vs CuPy).
  • Évolutivité : Le solveur a géré des grilles allant jusqu'à $120^3$ (environ 5 millions de DoF) efficacement.
  • Robustesse : La méthode est restée stable et a convergé de manière fluide même pour des matrices de permittivité mal conditionnées (valeurs propres allant de $10^{-1}$ à $10^{-5}$) violant les hypothèses standards, démontrant une robustesse supérieure aux approches traditionnelles.
• Nombre d'itérations : Le système anisotrope a nécessité légèrement plus d'itérations que le système isotrope, mais la méthode Cross-DoF a nécessité moins d'itérations que la méthode Triviale.

5. Importance

Cet article présente un solveur de pointe pour les cristaux photoniques dans les milieux anisotropes. Son importance réside dans :

• Combler le fossé : Il adapte avec succès le schéma de Yee efficace (généralement utilisé pour les milieux isotropes) à des cas anisotropes complexes avec une justification mathématique rigoureuse.
• Efficacité : En exploitant les opérations sans matrice et l'accélération GPU, il rend réalisables sur le plan computationnel les simulations 3D à grande échelle de structures photoniques complexes (comme les gyroïdes et les réseaux diamant).
• Fiabilité : Les garanties théoriques sur la propriété HPD et le coefficient de pénalité assurent que le solveur est robuste face aux instabilités numériques, en faisant un outil fiable pour la conception de dispositifs photoniques de nouvelle génération.

Le code source est disponible publiquement, facilitant la recherche et l'application ultérieures en photonique computationnelle.