Fast solvers for the high-order FEM simplicial de Rham complex: Extended edition

1. Problema e Motivação

O artigo aborda o desafio de resolver eficientemente os mapas de Riesz (e, consequentemente, os Laplacianos de Hodge) associados ao complexo de de Rham em malhas simpliciais (triângulos e tetraedros) utilizando elementos finitos de alta ordem ($p$).

• Desafio Computacional: Embora as discretizações de alta ordem ofereçam convergência rápida e alta intensidade aritmética, a implementação prática enfrenta barreiras de escalabilidade. Em uma célula $d$-dimensional, o número de graus de liberdade escala como $O(p^d)$. Com elementos padrão (tipo Lagrange), a montagem de matrizes densas e a resolução de sistemas lineares resultam em custos de $O(p^{3d})$ operações (flops) e $O(p^{2d})$ de armazenamento. Em 3D, isso significa $O(p^9)$ flops, tornando a abordagem ingênua proibitiva para altos valores de $p$.
• Limitação de Malhas Simpliciais: Em malhas tensoriais (hexaedros/quadriláteros), técnicas como fatoração de soma e bases tensoriais permitem reduzir o custo para $O(p^{d+1})$. No entanto, estender essas eficiências para malhas simpliciais (tetraedros) é significativamente mais difícil e, até então, não existiam solvers conhecidos com complexidade ótima para todo o complexo de de Rham nessas geometrias.

2. Metodologia

Os autores propõem uma nova abordagem baseada na construção de novos elementos finitos e estratégias de pré-condicionamento que exploram a estrutura algébrica do complexo de de Rham.

A. Novos Elementos Finitos e Graus de Liberdade

Em vez de usar bases padrão, os autores definem novos elementos baseados em uma ideia de Demkowicz et al., mas com uma escolha inteligente de graus de liberdade (GDLs) para promover ortogonalidade:

• Base no Simplexo de Referência: Os GDLs são definidos como momentos integrais sobre um simplexo de referência equilátero.
• Bases Ortogonais: Os GDLs utilizam uma base polinomial computada numericamente que é ortogonal em dois produtos internos diferentes: o produto interno $L^2$ e o produto interno associado ao operador diferencial ($H(dk)$).
• Decomposição Bubble: Os graus de liberdade são construídos resolvendo problemas de autovalor offline nas células de referência. Isso gera uma base onde as funções "interiores" (bubble functions) são fracamente acopladas com as funções de "interface" (bordas, faces, vértices).
• Propriedades da Matriz: Como resultado, as matrizes de massa e rigidez no interior das células tornam-se diagonais (ou quase diagonais), e o acoplamento interior-interior é nulo ou muito fraco.

B. Estratégias de Pré-condicionamento e Decomposição de Espaço

Para explorar a fraca acoplamento interior-interface, os autores combinam a nova base com métodos de decomposição de espaço (Schwarz):

1. Decoupling Interior-Interface: Eles separam os graus de liberdade em dois conjuntos:
   • Interior ($X_I$): Suportados apenas no interior de cada célula. Devido à diagonalização da base, estes podem ser resolvidos de forma barata usando pré-condicionadores do tipo Jacobi pontual (escalonamento diagonal).
   • Interface ($X_\Gamma$): Suportados nas fronteiras das células (vértices, arestas, faces).
2. Decomposições de Espaço:
   • PAFW (Pavarino-Arnold-Falk-Winther): Decomposição baseada em estrelas de vértices.
   • Hiptmair-Toselli (HT): Decomposição baseada em estrelas de entidades de dimensão $k$ e $k-1$ (mais eficiente para $k>0$).
   • Novidade: Os autores propõem uma variação onde apenas as funções do Tipo-I (relacionadas à imagem do operador diferencial) são usadas nas subproblemas de potencial, reduzindo drasticamente o tamanho dos subproblemas locais.
3. Algoritmo: O solver global utiliza um método de Schwarz aditivo/multiplicativo híbrido. As subproblemas de interface são resolvidos exatamente (ou com alta precisão), enquanto os graus de liberdade interiores são tratados com Jacobi pontual, evitando a necessidade de eliminação estática (static condensation) explícita e custosa.

3. Contribuições Chave

1. Complexidade Computacional Reduzida: O método proposto resolve os mapas de Riesz em 3D com complexidade de $O(p^6)$ flops e $O(p^4)$ armazenamento (na prática, $O(p^6)$ para a aplicação do operador e pré-condicionador), uma melhoria significativa em relação aos $O(p^9)$ das abordagens ingênuas. Embora não atinja o $O(p^4)$ ótimo das malhas tensoriais, é a primeira abordagem eficiente para malhas simpliciais de alta ordem.
2. Novos Elementos Finitos: A construção de bases que diagonalizam as submatrizes interiores de massa e rigidez no simplexo de referência, garantindo que o número de não-nulos por linha nas matrizes internas seja exatamente 1.
3. Estratégia de Pré-condicionamento Robusta: Demonstração teórica e numérica de que a separação interior-interface, combinada com decomposições de espaço (PAFW e HT), produz pré-condicionadores com contagens de iteração robustas em relação à ordem polinomial $p$ e ao parâmetro de malha $h$, desde que a razão $\beta h^2 / \alpha$ seja limitada.
4. Aplicação ao Laplaciano de Hodge: Extensão do método para pré-condicionar sistemas acoplados (Laplacianos de Hodge) usando pré-condicionadores de Lagrange aumentado, demonstrando eficácia em problemas com singularidades (como o canto de Fichera).

4. Resultados Numéricos

Os experimentos foram realizados usando a biblioteca Firedrake em malhas Freudenthal (refinadas) e no canto de Fichera.

• Mapas de Riesz ($H(\text{grad})$, $H(\text{curl})$, $H(\text{div})$):
  • Para $\alpha \ge 1$ (caso de interesse para pré-condicionadores de Lagrange aumentado), os solvers com separação interior-interface apresentam contagens de iteração quase idênticas aos solvers completos (sem separação), mas com custo computacional por iteração drasticamente menor.
  • A decomposição Hiptmair-Toselli baseada em Tipo-I (usando apenas funções de Tipo-I nos subespaços) mostrou-se a mais eficiente em termos de tamanho de subproblema, embora possa perder robustez em $p$ para casos onde $\alpha$ é muito pequeno ($\alpha \ll \beta$).
  • Para $H(\text{grad})$, a decomposição PAFW0 com separação interior-interface foi recomendada. Para $H(\text{curl})$, a decomposição HT com Tipo-I foi superior. Para $H(\text{div})$, uma variante baseada em estrelas de arestas (PAFW1) foi recomendada.
• Laplacianos de Hodge:
  • Ao aplicar os pré-condicionadores desenvolvidos em sistemas de Laplaciano de Hodge (problemas mistos), os autores observaram robustez em $p$ e convergência rápida.
  • O uso de $\gamma$ grande (Lagrange aumentado) reduziu significativamente o número de iterações (de ~50-60 para ~3-4 em alguns casos) sem custo adicional significativo no pré-condicionador.
  • A separação interior-interface manteve a eficiência, com contagens de iteração comparáveis às soluções completas.

5. Significado e Impacto

Este trabalho representa um avanço fundamental na computação científica de alta ordem para geometrias complexas (malhas não estruturadas de tetraedros).

• Viabilidade de Alta Ordem: Torna viável o uso de elementos finitos de alta ordem em malhas tetraédricas para problemas que exigem a estrutura do complexo de de Rham (como eletromagnetismo, mecânica dos fluidos e elasticidade), anteriormente limitados por custos computacionais proibitivos.
• Eficiência em Malhas Não Estruturadas: Preenche a lacuna entre a eficiência das malhas tensoriais e a flexibilidade geométrica das malhas simpliciais.
• Pré-condicionadores Modernos: Fornece a base para pré-condicionadores escaláveis para sistemas de equações diferenciais parciais complexos, permitindo a resolução de problemas em larga escala com alta precisão.
• Implementação Prática: A abordagem é implementável em pacotes modernos de elementos finitos (como Firedrake), facilitando a adoção por pesquisadores e engenheiros.

Em resumo, o artigo propõe uma solução elegante que combina a construção inteligente de bases de elementos finitos com métodos de domínio decomposto para superar a barreira de complexidade $O(p^9)$, permitindo solvers rápidos e robustos para o complexo de de Rham em 3D.