Geometric numerical discretization of electromagnetic quasineutral models

Enunciado do Problema
A simulação numérica de plasmas magnetizados é essencial para a energia de fusão, o clima espacial e a astrofísica, mas enfrenta desafios devido à vasta separação de escalas de comprimento e tempo. Embora modelos cinéticos baseados no sistema Vlasov-Maxwell sejam necessários para capturar distribuições não-maxwellianas e interações onda-partícula, resolver o comprimento de Debye e ondas eletromagnéticas de alta frequência (como ondas de Langmuir e de luz) é computacionalmente proibitivo para muitas aplicações. Modelos quase-neutros, que assumem neutralidade de carga ($\rho=0$) e filtram essas pequenas escalas, oferecem uma alternativa mais eficiente. No entanto, as abordagens numéricas existentes para sistemas quase-neutros frequentemente dependem de esquemas de preservação assintótica que requerem amortecimento numérico para suprimir modos de alta frequência conforme o comprimento de Debye se aproxima de zero, ou carecem da preservação da estrutura geométrica (como conservação de energia e restrições de divergência exatas) encontrada em métodos modernos de Particle-in-Cell (PIC) que preservam a estrutura.

Metodologia
Os autores estendem o framework de PIC eletromagnético geométrico, GEMPICX, para resolver as equações Vlasov-Maxwell totalmente cinéticas e quase-neutras. A metodologia é construída sobre os seguintes pilares:

• Framework Geométrico: O esquema utiliza diferenças finitas miméticas em grades duais (malhas primal e dual). As variáveis de campo são definidas em espaços discretos (nós, arestas, faces, volumes) que formam um complexo de de Rham discreto, garantindo a preservação exata de propriedades topológicas como a condição de divergência nula do campo magnético.
• Princípio Variacional: A discretização é derivada de um princípio de ação discreta baseado no Lagrangiano por Tronci & Camporeale (2015). Isso garante que as equações de movimento semi-discretas para partículas e campos sejam consistentes com a estrutura Hamiltoniana subjacente.
• Formulação Quase-neutra: Ao contrário das abordagens de preservação assintótica que resolvem o sistema Vlasov-Maxwell completo e tomam o limite $\epsilon_0 \to 0$, este método é derivado diretamente no limite quase-neutro. Consequentemente, o termo de corrente de deslocamento está ausente, e o campo elétrico não pode ser evoluído via uma derivada temporal. Em vez disso, o campo elétrico é obtido implicitamente em cada passo de tempo resolvendo uma equação curl-curl discretizada (Equação 2.25), que é não-singular e definida positiva.
• Correção de Divergência de Corrente: Para satisfazer a condição quase-neutra $\nabla \cdot \mathbf{J} = 0$ com precisão de máquina (o que não é automaticamente satisfeito pela deposição de partículas devido ao ruído estatístico), um multiplicador de Lagrange é introduzido. Isso se manifesta como um campo elétrico "corretor de divergência" ($\mathbf{E}^\star = -\nabla \psi^\star$) resolvido via uma equação do tipo Poisson (Equação 6.8) para ajustar as velocidades das partículas antes do passo de tempo principal.
• Passo de Tempo: O algoritmo emprega um esquema de divisão (split scheme) de segunda ordem de precisão. Ele envolve a correção da divergência da corrente, a resolução do potencial vetor $\mathbf{A}$ (usando o calibre de Coulomb $\nabla \cdot \mathbf{A} = 0$), a resolução do campo elétrico $\mathbf{E}$ e, finalmente, a atualização das posições e velocidades das partículas usando um Boris pusher.

Contribuições Principais

1. Extensão do GEMPICX: O artigo apresenta a primeira extensão do framework GEMPIC de diferenças finitas miméticas para o limite quase-neutro para elétrons e íons totalmente cinéticos.
2. Derivação Direta Quase-neutra: O esquema é derivado diretamente no limite quase-neutro em vez de depender da consistência assintótica de um solver completo. Isso garante a ausência de ondas eletromagnéticas e de Langmuir espúrias na solução numérica sem depender de amortecimento numérico.
3. Preservação de Estrutura: A formulação semi-discreta prova matematicamente a conservação da soma das energias magnética e cinética. Ela também mantém o requisito de divergência zero para a densidade de corrente com precisão de máquina através do uso de um multiplicador de Lagrange.
4. Solver de Campo Elétrico Implícito: Os autores derivam um sistema linear robusto (Equação 5.12) para obter o campo elétrico implicitamente, contornando a necessidade de potenciais elétricos na atualização do campo e evitando os problemas de singularidade associados ao operador curl-curl padrão na ausência de corrente de deslocamento.
5. Escalabilidade: A implementação é integrada ao framework AMReX com paralelização MPI, demonstrando capacidades de escalonamento forte e fraco em recursos de HPC.

Resultados
O framework foi validado através de uma série de testes numéricos com condições de contorno periódicas:

• Relações de Dispersão: Simulações de plasma de espécie única (apenas elétrons) e de duas espécies (elétrons-íons) reproduziram com sucesso as relações de dispersão teóricas para ondas propagando-se paralelamente e perpendicularmente ao campo magnético de fundo.
  • Propagação Paralela: O esquema capturou corretamente a onda de ciclotron de elétrons (ECW) e a onda de ciclotron de íons (ICW), bem como as estruturas cônicas de modos Alfvén-ciclotron de ordem superior amortecidos. Como esperado, as ondas de Langmuir e os modos eletromagnéticos L/R estavam ausentes.
  • Propagação Perpendicular: As simulações reproduziram ondas de Bernstein e ressonâncias harmônicas de ciclotron. Em casos de duas espécies, foram observadas ondas de Bernstein de íons e a onda híbrida. Os modos O e X estavam corretamente ausentes devido à suposição quase-neutra.
• Amortecimento de Onda: Um teste envolvendo uma perturbação senoidal única demonstrou a capacidade do código em capturar a taxa de amortecimento prevista teoricamente de uma onda de ciclotron, validando a precisão do esquema no tratamento de processos cinéticos de ordem superior.
• Eficiência Computacional: Uma comparação entre a resolução do sistema curl-curl completo para o campo elétrico (Equação 5.12) e uma aproximação diagonal mais barata (Equação 5.13) mostrou que a aproximação reduziu o custo computacional em aproximadamente 25% com perda negligível de precisão para os casos testados.
• Escalabilidade: Estudos de escalonamento forte e fraco em uma arquitetura de memória distribuída mostraram um aumento de velocidade (speedup) quase ideal até 32 tarefas MPI e mantiveram alta eficiência até 12 mais 128 tarefas, indicando que o código é bem adequado para simulações de larga escala.

Significância e Alegações
Os autores alegam que este trabalho consegue preencher a lacuna entre métodos de preservação de estrutura geométrica e modelagem de plasma quase-neutro. Ao derivar o esquema diretamente no limite quase-neutro, o método filtra inerentemente modos de alta frequência, oferecendo uma alternativa mais eficiente e fisicamente consistente aos esquemas de preservação assintótica que dependem de amortecimento. A preservação da estrutura de de Rham discreta e da conservação de energia (no nível semi-discreto) garante estabilidade de longo prazo e fidelidade física.

O artigo nota modestamente limitações e áreas para trabalhos futuros:

• Condições de Contorno: A implementação atual depende de condições de contorno periódicas para impor facilmente o calibre de Coulomb e inverter o operador de Laplace. Estender isso para contornos não periódicos (por exemplo, para dispositivos de fusão) requer análises de compatibilidade e restrições adicionais.
• Ordem de Precisão: O trabalho utiliza operadores de Hodge de segunda ordem. Embora operadores de ordem superior sejam possíveis dentro do framework mimético, seu benefício específico para simulações quase-neutras (que carecem de modos de alta frequência) permanece uma questão aberta.
• Tratamento de Choques: O método assume campos suaves. Embora os métodos PIC geralmente lidem com choques através de distribuições de partículas, o comportamento desta discretização geométrica específica em regimes com gradientes fortes ou choques requer investigação adicional.
• Discretização Temporal: A conservação de energia é exata apenas no nível semi-discreto; após a discretização temporal, a conservação é precisa apenas até a ordem do integrador de tempo. Uma discretização temporal totalmente variacional é sugerida como uma extensão futura para alcançar a conservação de energia discreta exata.