Geometric numerical discretization of electromagnetic quasineutral models

Énoncé du problème
La simulation numérique des plasmas magnétisés est essentielle pour la fusion nucléaire, la météo spatiale et l'astrophysique, mais elle se heurte à des défis dus à la vaste séparation des échelles de longueur et de temps. Bien que les modèles cinétiques basés sur le système Vlasov-Maxwell soient nécessaires pour capturer les distributions non maxwelliennes et les interactions onde-particule, la résolution de la longueur de Debye et des ondes électromagnétiques à haute fréquence (telles que les ondes de Langmuir et les ondes lumineuses) est informatiquement prohibitive pour de nombreuses applications. Les modèles quasi-neutres, qui supposent la neutralité de charge ($\rho=0$) et filtrent ces petites échelles, offrent une alternative plus efficace. Cependant, les approches numériques existantes pour les systèmes quasi-neutres reposent souvent sur des schémas préservant l'asymptotique qui nécessitent un amortissement numérique pour supprimer les modes à haute fréquence à mesure que la longueur de Debye tend vers zéro, ou manquent de la préservation de la structure géométrique (telle que la conservation de l'énergie et les contraintes de divergence exacte) propre aux méthodes modernes de Particle-in-Cell (PIC) préservant la structure.

Méthodologie
Les auteurs étendent le cadre GEMPICX, un cadre de PIC électromagnétique géométrique, pour résoudre les équations de Vlasov-Maxwell pleinement cinétiques et quasi-neutres. La méthodologie repose sur les piliers suivants :

• Cadre géométrique : Le schéma utilise des différences finies mimétiques sur des grilles duales (maillages primaux et duaux). Les variables de champ sont définies sur des espaces discrets (nœuds, arêtes, faces, volumes) qui forment un complexe de de Rham discret, garantissant la préservation exacte des propriétés topologiques comme la condition de divergence nulle du champ magnétique.
• Principe variationnel : La discrétisation est dérivée d'un principe d'action discret basé sur le Lagrangien par Tronci & Camporeale (2015). Cela garantit que les équations du mouvement semi-discrètes pour les particules et les champs sont cohérentes avec la structure hamiltonienne sous-jacente.
• Formulation quasi-neutre : Contrairement aux approches préservant l'asymptotique qui résolvent le système complet Vlasov-Maxwell et prennent la limite $\epsilon_0 \to 0$, cette méthode est dérivée directement dans la limite quasi-neutre. Par conséquent, le terme de courant de déplacement est absent et le champ électrique ne peut pas être évolué via une dérivée temporelle. Au lieu de cela, le champ électrique est obtenu implicitement à chaque pas de temps en résolvant une équation curl-curl discrétisée (Équation 2.25), qui est non singulière et définie positive.
• Correction de la divergence du courant : Pour satisfaire la condition de quasi-neutralité $\nabla \cdot \mathbf{J} = 0$ à la précision machine (ce qui n'est pas automatiquement satisfait par le dépôt de particules à cause du bruit statistique), un multiplicateur de Lagrange est introduit. Cela se manifeste par un champ électrique « correcteur de divergence » ($\mathbf{E}^\star = -\nabla \psi^\star$) résolu via une équation de type Poisson (Équation 6.8) pour ajuster les vitesses des particules avant le pas de temps principal.
• Pas de temps : L'algorithme emploie un schéma fractionné (split scheme) d'ordre deux. Il implique la correction de la divergence du courant, la résolution du potentiel vecteur $\mathbf{A}$ (en utilisant la jauge de Coulomb $\nabla \cdot \mathbf{A} = 0$), la résolution du champ électrique $\mathbf{E}$, et enfin la mise à jour des positions et des vitesses des particules à l'aide d'un pousseur de Boris (Boris pusher).

Contributions clés

1. Extension de GEMPICX : L'article présente la première extension du cadre de différences finies mimétiques GEMPIC à la limite quasi-neutre pour les électrons et les ions pleinement cinétiques.
2. Dérivation directe quasi-neutre : Le schéma est dérivé directement dans la limite quasi-neutre plutôt que de compter sur la cohérence asymptotique d'un solveur complet. Cela garantit l'absence d'ondes électromagnétiques et de Langmuir parasites dans la solution numérique sans dépendre d'un amortissement numérique.
3. Préservation de la structure : La formulation semi-discrète prouve mathématiquement la conservation de la somme des énergies magnétique et cinétique. Elle maintient également l'exigence de divergence nulle de la densité de courant à la précision machine grâce à l'utilisation d'un multiplicateur de Lagrange.
4. Solveur de champ électrique implicite : Les auteurs dérivent un système linéaire robuste (Équation 5.12) pour obtenir le champ électrique implicitement, évitant ainsi le recours aux potentiels électriques pour la mise à jour du champ et contournant les problèmes de singularité associés à l'opérateur curl-curl standard en l'absence de courant de déplacement.
5. Scalabilité : L'implémentation est intégrée au framework AMReX avec une parallélisation MPI, démontrant des capacités de mise à l'échelle forte et faible.

Résultats
Le cadre a été validé par une série de tests numériques avec des conditions aux limites périodiques :

• Relations de dispersion : Les simulations de plasma à espèce unique (électron seul) et à double espèce (électron-ion) ont reproduit avec succès les relations de dispersion théoriques pour les ondes se propageant parallèlement et perpendiculairement au champ magnétique de fond.
  • Propagation parallèle : Le schéma a correctement capturé l'onde cyclotronique électronique (ECW) et l'onde cyclotronique ionique (ICW), ainsi que les structures coniques des modes Alfvén-cyclotron d'ordre supérieur amortis. Comme prévu, les ondes de Langmuir et les modes électromagnétiques L/R étaient absents.
  • Propagation perpendiculaire : Les simulations ont reproduit les ondes de Bernstein et les résonances harmoniques cyclotroniques. Dans les cas à double espèce, les ondes de Bernstein ioniques et l'onde hybride ont été observées. Les modes O et X étaient correctement absents en raison de l'hypothèse de quasi-neutralité.
• Amortissement des ondes : Un test impliquant une perturbation sinusoïdale unique a démontré la capacité du code à capturer le taux d'amortissement théoriquement prédit d'une onde cyclotronique, validant la précision du schéma dans la gestion des processus cinétiques d'ordre supérieur.
• Efficacité de calcul : Une comparaison entre la résolution du système curl-curl complet pour le champ électrique (Équation 5.12) et une approximation diagonale moins coûteuse (Équation 5.13) a montré que l'approximation réduisait le coût de calcul d'environ 25 % avec une perte de précision négligeable pour les cas testés.
• Scalabilité : Les études de mise à l'échelle forte et faible sur une architecture à mémoire distribuée ont montré une accélération quasi idéale jusqu'à 32 tâches MPI et ont maintenu une haute efficacité jusqu'à 128 tâches, indiquant que le code est bien adapté aux simulations à grande échelle.

Signification et revendications
Les auteurs affirment que ce travail comble avec succès le fossé entre les méthodes géométriques préservant la structure et la modélisation de plasma quasi-neutre. En dérivant le schéma directement dans la limite quasi-neutre, la méthode filtre intrinsèquement les modes à haute fréquence, offrant une alternative plus efficace et physiquement cohérente aux schémas préservant l'asymptotique qui reposent sur l'amortissement. La préservation de la structure de de Rham discrète et de la conservation de l'énergie (au niveau semi-discret) assure la stabilité à long terme et la fidélité physique.

Le papier note modestement les limites et les domaines de travaux futurs :

• Conditions aux limites : L'implémentation actuelle repose sur des conditions aux limites périodiques pour imposer facilement la jauge de Coulomb et inverser l'opérateur de Laplace. L'extension à des limites non périodiques (par exemple, pour les dispositifs de fusion) nécessite des contraintes supplémentaires et une analyse de compatibilité.
• Ordre de précision : Le travail utilise des opérateurs de Hodge d'ordre deux. Bien que des opérateurs d'ordre supérieur soient possibles dans le cadre mimétique, leur bénéfice spécifique pour les simulations quasi-neutres (qui sont dépourvues de modes à haute fréquence) reste une question ouverte.
• Gestion des chocs : La méthode suppose des champs lisses. Bien que les méthodes PIC gèrent généralement les chocs via les distributions de particules, le comportement de cette discrétisation géométrique spécifique dans des régimes avec de forts gradients ou des chocs nécessite des investigations supplémentaires.
• Discrétisation temporelle : La conservation de l'énergie est exacte uniquement au niveau semi-discret ; après la discrétisation temporelle, la conservation n'est précise qu'à l'ordre de l'intégrateur temporel. Une discrétisation temporelle entièrement variationnelle est suggérée comme une extension future pour atteindre une conservation d'énergie discrète exacte.