A dimension-oblivious domain decomposition method based on space-filling curves

1. Formulação do Problema

O artigo aborda o desafio de resolver equações diferenciais parciais (EDPs) elípticas de alta dimensão, especificamente a equação de Poisson, em sistemas modernos de computação exascale. As principais dificuldades identificadas são:

• Maldição da Dimensionalidade: Métodos de discretização padrão (malhas completas) tornam-se inviáveis à medida que a dimensão $d$ aumenta, devido ao crescimento exponencial dos graus de liberdade (GDLs).
• Malhas Anisotrópicas: A "Técnica de Combinação de Malhas Esparsas" é utilizada para mitigar a maldição da dimensionalidade. Este método decompõe o problema em muitos subproblemas independentes em malhas anisotrópicas (onde os tamanhos de malha $h_j$ variam por dimensão). Métodos tradicionais de decomposição de domínio (DD) geométrica lutam aqui porque dependem de informações geométricas de dimensão fixa, o que é difícil de definir para malhas anisotrópicas arbitrárias.
• Tolerância a Falhas e Escalabilidade: Sistemas exascale futuros exigem algoritmos que não sejam apenas escaláveis para milhões de núcleos, mas também tolerantes a falhas. Atingir tolerância a falhas em solvers numéricos geralmente requer redundância de dados, o que implica grandes sobreposições entre subdomínios em um método DD. Métodos DD convencionais minimizam a sobreposição para reduzir custos de comunicação, criando um conflito com a necessidade de redundância.
• Ignorância de Dimensão: O solver deve operar algebricamente nas equações matriciais montadas, sem depender de coordenadas geométricas, permitindo lidar com dimensões arbitrárias e malhas anisotrópicas de forma uniforme.

2. Metodologia

Os autores propõem um solver de decomposição de domínio algébrico de dois níveis baseado em Curvas de Preenchimento de Espaço (SFCs).

A. Particionamento via Curvas de Preenchimento de Espaço

Em vez de particionamento geométrico, o método mapeia os pontos da malha $d$-dimensional para um espaço de parâmetro unidimensional $[0, 1]$ usando uma curva de preenchimento de espaço discreta (especificamente a curva de Hilbert).

• Mapeamento: Os pontos da malha são ordenados de acordo com a SFC.
• Particionamento Disjunto: A sequência ordenada de $N$ pontos é dividida em $P$ subconjuntos contíguos de tamanho aproximadamente igual ($N/P$). Isso garante uma carga balanceada independentemente da dimensão $d$ ou da anisotropia da malha.
• Construção de Sobreposição: Para alcançar redundância de dados para tolerância a falhas, os subconjuntos disjuntos são ampliados. Em vez de adicionar camadas de fronteira geométricas, o método adiciona conjuntos de índices vizinhos da ordenação 1D da SFC. Um parâmetro de sobreposição $\gamma$ controla quantos subdomínios vizinhos são incluídos.
  • Para $\gamma = 0.5$, cada ponto da malha é coberto exatamente por $2\gamma + 1 = 2$ subdomínios.
  • Para múltiplos inteiros de $0.5$, a sobreposição é uniforme, garantindo simetria nos operadores resultantes.

B. Solver Algébrico de Dois Níveis

O solver é construído puramente de forma algébrica, evitando malhas grosseiras geométricas:

1. Nível Fino (Subdomínios Locais): Cada processador resolve um problema local em seu subdomínio sobreposto $\Omega_i$.
2. Nível Grosseiro (Correção Global): Um espaço grosseiro é definido algebricamente. Em vez de uma malha grosseira geométrica, o espaço grosseiro é construído atribuindo $q$ graus de liberdade por subdomínio com base na ordenação da SFC. A matriz grosseira $A_0$ é formada via abordagem de Galerkin ($A_0 = R_0 A R_0^T$).
3. Redundância (Técnica Bank-Holst): O problema grosseiro é resolvido redundantemente em cada processador, juntamente com o subproblema local. Isso evita gargalos de comunicação associados a uma única solução global grosseira e suporta tolerância a falhas.
4. Ponderação: Para lidar com a contagem excessiva de incógnitas em regiões sobrepostas, aplica-se uma partição da unidade. Os autores mostram que, para escolhas específicas de $\gamma$ (múltiplos de $0.5$), as matrizes de ponderação $D_i$ tornam-se múltiplos escalares da identidade, preservando a simetria do operador. Isso permite o uso do método do Gradiente Conjugado (CG).

C. Esquemas Iterativos

O método é implementado como:

• Uma iteração de Richardson amortecida (Schwarz Aditivo).
• Um método de Gradiente Conjugado Precondicionado (PCG).
• Ambos utilizam o operador de dois níveis (Soluções locais + Correção grosseira).

3. Contribuições Principais

• Solver Ignorante de Dimensão: O primeiro solver DD que opera diretamente em matrizes montadas usando SFCs, tornando-o aplicável a qualquer dimensão $d$ e qualquer estrutura de malha anisotrópica sem suposições geométricas.
• Grande Sobreposição para Tolerância a Falhas: O artigo projeta intencionalmente o método DD com grandes sobreposições (controladas por $\gamma$) para fornecer a redundância de dados necessária para computação tolerante a falhas em sistemas exascale, um desvio do DD tradicional que minimiza a sobreposição.
• Malha Grosseira Algébrica: Uma construção puramente algébrica do espaço grosseiro que evita a complexidade de definir malhas grosseiras geométricas para problemas anisotrópicos ou de alta dimensão.
• Preservação de Simetria: Prova teórica (Lema 3.1) de que escolher $\gamma$ como um múltiplo de $0.5$ resulta em operadores simétricos, permitindo o uso de solvers simétricos eficientes como o PCG, mesmo com grandes sobreposições.

4. Resultados Numéricos

Os autores conduziram extensos experimentos numéricos em dimensões $d=1$ a $6$ usando até um milhão de processadores.

• Comportamento de Convergência:

  • O solver exibe escalabilidade fraca ótima: O número de iterações necessárias para convergir permanece constante à medida que o número de processadores $P$ aumenta (mantendo o tamanho do problema por processador fixo).
  • Independência de Dimensão: A taxa de convergência é amplamente independente da dimensão $d$. Curiosamente, o caso 1D é o mais difícil (convergência mais lenta) devido à distância relativa entre as escalas fina e grosseira; dimensões mais altas convergem na verdade mais rápido ou de forma similar.
  • PCG vs. Richardson: O método de Gradiente Conjugado Precondicionado supera significativamente a iteração de Richardson, exigindo aproximadamente metade das iterações.
  • Parâmetro de Sobreposição: Uma sobreposição de $\gamma = 0.5$ fornece o melhor equilíbrio entre velocidade de convergência e redundância. Maior $\gamma$ melhora a tolerância a falhas, mas aumenta ligeiramente as contagens de iteração (a menos que ponderado corretamente).

• Escalabilidade na Técnica de Combinação:

  • O método foi integrado à Técnica de Combinação de Malhas Esparsas.
  • Paralelismo Massivo: A abordagem utilizou com sucesso mais de 1 milhão de processadores para um problema de 6 dimensões.
  • Balanceamento de Carga: O particionamento SFC garantiu que cada processador lidasse com aproximadamente o mesmo número de GDLs ($2^8$), resultando em tempos de execução altamente uniformes em todos os subproblemas.
  • Eficiência: O tempo de execução por subproblema permaneceu constante à medida que o tamanho total do sistema crescia, demonstrando verdadeira escalabilidade fraca.

5. Significado

Este trabalho fornece um bloco de construção crítico para a computação exascale em EDPs de alta dimensão.

• Escalabilidade: Resolve a "maldição da dimensionalidade" não apenas usando malhas esparsas, mas permitindo que a técnica de combinação de malhas esparsas seja executada em milhões de núcleos de forma eficiente.
• Tolerância a Falhas: Ao deslocar o foco de minimizar a comunicação (pequena sobreposição) para maximizar a redundância (grande sobreposição), o método abre caminho para solvers que podem sobreviver a falhas de hardware sem reiniciar toda a simulação.
• Generalidade: A natureza algébrica e ignorante de dimensão do solver significa que ele pode ser aplicado a uma ampla gama de problemas (elípticos, parabólicos) e discretizações (diferenças finitas, elementos finitos) sem necessidade de ajuste geométrico específico para dimensão.

Em resumo, o artigo apresenta um solver linear robusto, escalável e tolerante a falhas que aproveita curvas de preenchimento de espaço para desacoplar o processo de decomposição de domínio da dimensão geométrica do problema, permitindo a solução eficiente de EDPs de alta dimensão em futuros supercomputadores.