Geometric numerical discretization of electromagnetic quasineutral models

Problematica
La simulazione numerica di plasmi magnetizzati è essenziale per la fusione nucleare, il meteo spaziale e l'astrofisica, tuttavia affronta sfide dovute alla vasta separazione delle scale di lunghezza e di tempo. Mentre i modelli cinetici basati sul sistema Vlasov-Maxwell sono necessari per catturare distribuzioni non-Maxwelliane e interazioni onda-particella, la risoluzione della lunghezza di Debye e delle onde elettromagnetiche ad alta frequenza (come le onde di Langmuir e le onde luminose) è computazionalmente proibitiva per molte applicazioni. I modelli quasineutrali, che assumono la neutralità della carica ($\rho=0$) e filtrano queste piccole scale, offrono un'alternativa più efficiente. Tuttavia, gli approcci numerici esistenti per sistemi quasineutrali spesso si affidano a schemi asintoticamente preservanti (asymptotic-preserving) che richiedono uno smorzamento numerico per sopprimere i modi ad alta frequenza man mano che la lunghezza di Debye tende a zero, oppure mancano della preservazione della struttura geometrica (come la conservazione dell'energia e i vincoli esatti di divergenza) tipica dei moderni metodi Particle-in-Cell (PIC) che preservano la struttura.

Metodologia
Gli autori estendono il framework geometrico elettrodinamico GEMPICX per risolvere le equazioni Vlasov-Maxwell completamente cinetiche e quasineutrali. La metodologia si basa sui seguenti pilastri:

• Framework Geometrico: Lo schema utilizza differenze finite mimetiche su griglie duali (mesh primale e duale). Le variabili di campo sono definite su spazi discreti (nodi, spigoli, facce, volumi) che formano un complesso de Rham discreto, garantendo la preservazione esatta delle proprietà topologiche come la condizione di divergenza nulla del campo magnetico.
• Principio Variazionale: La discretizzazione è derivata da un principio d'azione discreto basato sul Lagrangiano di Tronci & Camporeale (2015). Ciò assicura che le equazioni del moto semi-discrete per particelle e campi siano coerenti con la sottostante struttura Hamiltoniana.
• Formulazione Quasineutrale: A differenza degli approcci asymptotic-preserving che risolvono il sistema completo Vlasov-Maxwell e prendono il limite $\epsilon_0 \to 0$, questo metodo è derivato direttamente nel limite quasineutrale. Di conseguenza, il termine di corrente di spostamento è assente e il campo elettrico non può essere evoluto tramite una derivata temporale. Invece, il campo elettrico viene ottenuto implicitamente ad ogni passo temporale risolvendo un'equazione curl-curl discretizzata (Equazione 2.25), che è non singolare e definita positiva.
• Correzione della Divergenza di Corrente: Per soddisfare la condizione di quasineutralità $\nabla \cdot \mathbf{J} = 0$ con precisione di macchina (il che non è automaticamente soddisfatto dalla deposizione delle particelle a causa del rumore statistico), viene introdotto un moltiplicatore di Lagrange. Questo si manifesta come un campo elettrico "correttore di divergenza" ($\mathbf{E}^\star = -\nabla \psi^\star$) risolto tramite un'equazione di tipo Poisson (Equazione 6.8) per regolare le velocità delle particelle prima del passo temporale principale.
• Integrazione Temporale (Time-Stepping): L'algoritmo impiega uno schema split di secondo ordine. Esso prevede la correzione della divergenza di corrente, la risoluzione del potenziale vettore $\mathbf{A}$ (utilizzando il gauge di Coulomb $\nabla \cdot \mathbf{A} = 0$), la risoluzione del campo elettrico $\mathbf{E}$ e, infine, l'aggiornamento delle posizioni e delle velocità delle particelle utilizzando un pusher di Boris.

Contributi Chiave

1. Estensione di GEMPICX: Il lavoro presenta la prima estensione del framework GEMPIC a differenze finite mimetiche al limite quasineutrale per elettroni e ioni completamente cinetici.
2. Derivazione Diretta Quasineutrale: Lo schema è derivato direttamente nel limite quasineutrale anziché affidarsi alla coerenza asintotica di un solver completo. Ciò garantisce l'assenza di onde elettromagnetiche e di Langmuir spurie nella soluzione numerica senza ricorrere allo smorzamento numerico.
3. Preservazione della Struttura: La formulazione semi-discreta dimostra matematicamente la conservazione della somma delle energie magnetica e cinetica. Essa mantiene inoltre il requisito di divergenza nulla della densità di corrente con precisione di macchina attraverso l'uso di un moltiplicatore di Lagrange.
4. Solver Implicito del Campo Elettrico: Gli autori derivano un sistema lineare robusto (Equazione 5.12) per ottenere il campo elettrico implicitamente, evitando la necessità di potenziali elettrici nell'aggiornamento del campo e scongiurando i problemi di singolarità associati all'operatore curl-curl standard in assenza di corrente di spostamento.
5. Scalabilità: L'implementazione è integrata nel framework AMReX con parallelizzazione MPI, dimostrando capacità di scalabilità forte e debole su risorse HPC.

Risultati
Il framework è stato validato attraverso una serie di test numerici con condizioni al contorno periodiche:

• Relazioni di Dispersione: Le simulazioni di plasma a specie singola (solo elettroni) e a doppia specie (elettroni-ioni) hanno riprodotto con successo le relazioni di dispersione teoriche per onde che si propagano parallelamente e perpendicolarmente al campo magnetico di fondo.
  • Propagazione Parallela: Lo schema ha catturato correttamente l'onda ciclotronica elettronica (ECW) e l'onda ciclotronica ionica (ICW), nonché le strutture coniche dei modi Alfvén-ciclotronici di ordine superiore smorzati. Come previsto, le onde di Langmuir e le onde elettromagnetiche L/R erano assenti.
  • Propagazione Perpendicolare: Le simulazioni hanno riprodotto le onde di Bernstein e le risonanze armoniche ciclotroniche. Nei casi a doppia specie, sono state osservate le onde di Bernstein ioniche e l'onda ibrida. Le modalità O-mode e X-mode erano correttamente assenti a causa dell'assunzione di quasineutralità.
• Smorzamento delle Onde: Un test che coinvolgeva una singola perturbazione sinusoidale ha dimostrato la capacità del codice di catturare il tasso di smorzamento teoricamente previsto di un'onda ciclotronica, validando l'accuratezza dello schema nel gestire processi cinetici di ordine superiore.
• Efficienza Computazionale: Un confronto tra la risoluzione del sistema curl-curl completo per il campo elettrico (Equazione 5.12) e un'approssimazione diagonale più economica (Equazione 5.13) ha mostrato che l'approssimazione riduce il costo computazionale di circa il 25% con una perdita di accuratezza trascurabile per i casi testati.
• Scalabilità: Studi di scalabilità forte e debole su un'architettura a memoria distribuita hanno mostrato un incremento di velocità (speedup) quasi ideale fino a 32 task MPI e hanno mantenuto un'alta efficienza fino a 128 task, indicando che il codice è ben adatto per simulazioni su larga scala.

Significatività e Rivendicazioni
Gli autori affermano che questo lavoro riesce a colmare il divario tra i metodi geometrici che preservano la struttura e la modellazione di plasmi quasineutrali. Derivando lo schema direttamente nel limite quasineutrale, il metodo filtra intrinsecamente i modi ad alta frequenza, offrendo un'alternativa più efficiente e fisicamente coerente rispetto agli schemi asymptotic-preserving che si affidano allo smorzamento. La preservazione della struttura de Rham discreta e della conservazione dell'energia (al livello semi-discreto) garantisce stabilità a lungo termine e fedeltà fisica.

Il documento nota con modestia i limiti e le aree di lavoro futuro:

• Condizioni al Contorno: L'attuale implementazione si affida a condizioni al contorno periodiche per imporre facilmente il gauge di Coulomb e invertire l'operatore di Laplace. Estendere questo a frontiere non periodiche (ad esempio, per dispositivi di fusione) richiede vincoli aggiuntivi e un'analisi di compatibilità.
• Ordine di Accuratezza: Il lavoro utilizza operatori di Hodge del secondo ordine. Sebbene siano possibili operatori di ordine superiore all'interno del framework mimetico, il loro specifico beneficio per le simulazioni quasineutrali (che mancano di modi ad alta frequenza) rimane una questione aperta.
• Gestione degli Shock: Il metodo assume campi regolari. Sebbene i metodi PIC gestiscano generalmente gli shock attraverso le distribuzioni delle particelle, il comportamento di questa specifica discretizzazione geometrica in regimi con forti gradienti o shock richiede ulteriori indagini.
• Discretizzazione Temporale: La conservazione dell'energia è esatta solo al livello semi-discreto; dopo la discretizzazione temporale, la conservazione è accurata solo fino all'ordine dell'integratore temporale. Una discretizzazione temporale completamente variazionale è suggerita come estensione futura per ottenere una conservazione dell'energia discreta esatta.