Cálculo da estrutura de escarsidade para matrizes de elementos finitos

13

Pergunta: Quais métodos estão disponíveis para calcular com precisão e eficiência a estrutura de esparsidade de uma matriz de elementos finitos?

Info: Estou trabalhando em um solucionador de equações de pressão de Poisson, usando o método de Galerkin com base em Lagrange quadrático, escrito em C, e usando o PETSc para armazenamento de matriz esparsa e rotinas KSP. Para usar o PETSc com eficiência, preciso pré-alocar memória para a matriz de rigidez global.

Atualmente, estou fazendo uma montagem simulada para estimar o número de não-zeros por linha da seguinte forma (pseudocódigo)

int nnz[global_dim]
for E=1 to NUM_ELTS
  for i=1 to 6
    gi = global index of i 
    if node gi is free
      for j=1 to 6
        gj = global index of j
        if node gj is free 
          nnz[i]++

No entanto, isso superestima o nnz porque algumas interações nó-nó podem ocorrer em vários elementos.

Pensei em tentar rastrear quais interações j, eu encontrei, mas não tenho certeza de como fazer isso sem usar muita memória. Eu também poderia fazer um loop sobre os nós e encontrar o suporte da função base centralizada nesse nó, mas então eu teria que pesquisar todos os elementos para cada nó, o que parece ineficiente.

Encontrei esta pergunta recente, que continha algumas informações úteis, especialmente de Stefano M, que escreveu

meu conselho é implementá-lo em python ou C, aplicando alguns conceitos teóricos dos grafos, ou seja, considere os elementos da matriz como arestas em um gráfico e calcule a estrutura de esparsidade da matriz de adjacência. Lista de listas ou dicionário de chaves são escolhas comuns.

Estou procurando mais detalhes e recursos sobre isso. Admito que não conheço muita teoria dos grafos e não estou familiarizado com todos os truques de CS que podem ser úteis (estou abordando isso do lado matemático).

Obrigado!

matrix finite-element petsc performance John Edwardson
fonte

5

Sua idéia de acompanhar quais interações i, j que você encontrou pode funcionar, acho que esse é o "truque do CS" ao qual você e Stefano M estão se referindo. Isso equivale a construir sua matriz esparsa no formato de lista de listas .

Não tenho certeza de quanto CS você tem, portanto peço desculpas se isso já é conhecido: em uma estrutura de dados de lista vinculada , cada entrada armazena um ponteiro para a entrada depois dela e a entrada anterior. É barato adicionar e excluir entradas de, mas não é tão simples encontrar itens nele - talvez você precise examinar todos eles.

Portanto, para cada nó i, você armazena uma lista vinculada. Então você itera através de todos os elementos; se você encontrar dois nós iej estão conectados, procure na lista vinculada de i. Se j ainda não estiver lá, adicione-o à lista e da mesma forma adicione i à lista de j. É mais fácil se você adicioná-los em ordem.

Depois de preencher sua lista de listas, agora você sabe o número de entradas diferentes de zero em cada linha da matriz: é o comprimento da lista desse nó. Essas informações são exatamente o que você precisa para pré-alocar uma matriz esparsa na estrutura de dados da matriz do PETSc. Depois, você pode liberar sua lista de listas porque não precisa mais dela.

No entanto, essa abordagem supõe que tudo o que você tem é a lista de quais nós cada elemento contém.

Alguns pacotes de geração de malha - Triângulo por exemplo - podem gerar não apenas uma lista de elementos e quais nós eles contêm, mas também uma lista de todas as arestas da triangulação. Nesse caso, você não corre o risco de superestimar o número de entradas diferentes de zero: para elementos lineares por partes, cada aresta fornece exatamente 2 entradas da matriz de rigidez. Você está usando quadrático por partes, de modo que cada aresta conta com 4 entradas, mas você entendeu. Nesse caso, você pode encontrar o número de entradas diferentes de zero por linha com uma passagem na lista de arestas usando uma matriz comum.

Com essa abordagem, é necessário ler um arquivo grande extra do disco rígido, que pode ser mais lento do que usar a lista de elementos se o cálculo real não for tão grande. No entanto, acho que é mais simples.

Daniel Shapero
fonte

Obrigado. Eu tenho uma lista de borda disponível, então provavelmente usarei seu segundo método por enquanto, mas posso voltar e tentar o primeiro método, apenas para sujar minhas mãos com listas vinculadas e outras coisas (obrigado pela introdução ... já única tomado uma classe básica CS, e enquanto eu faço talento especial para a programação, eu não sei tanto quanto eu deveria sobre estruturas de dados e algoritmos)

John Edwardson

Feliz por ajudar! Peguei muito do meu conhecimento de CS a partir disso: books.google.com/books?isbn=0262032937 - pelo amor de Deus, leia sobre análise amortizada. Programar sua própria lista vinculada ou estrutura de dados da árvore de pesquisa binária em C vale a pena.

Daniel Shapero

5

Se você especificar sua malha como um DMPlex e seu layout de dados como um PetscSection, o DMCreateMatrix () fornecerá automaticamente a matriz pré-alocada corretamente. Aqui estão exemplos do PETSc para o problema de Poisson e o problema de Stokes .

Matt Knepley
fonte

2

Voto a favor

Eu pessoalmente não conheço nenhuma maneira barata de fazer isso, então simplesmente superestimo o número, ou seja, uso um valor razoavelmente grande para todas as linhas.

Por exemplo, para uma malha perfeitamente estruturada feita de elementos hexagonais lineares de 8 nós, os nnzs por linha nos blocos diagonais e diagonais externos são dof * 27. Para a maioria das malhas sextavadas geradas automaticamente não estruturadas, o número raramente excede dof * 54. Para tets lineares, nunca tive a necessidade de ir além de dof * 30. Para algumas malhas com elementos de formato muito baixo / com formato muito baixo, pode ser necessário usar valores um pouco maiores.

A penalidade é que o consumo de memória local (na classificação) está entre 2x-5x, portanto, talvez você precise usar mais nós de computação no cluster do que o habitual.

Btw eu tentei usar listas pesquisáveis, mas o tempo necessário para determinar a estrutura de escarsidade foi mais do que a montagem / resolução. Mas minha implementação foi muito simples e não utilizou informações sobre arestas.

A outra opção é usar rotinas como DMMeshCreateExodus, como mostrado neste exemplo.

stali
fonte

0

Você está procurando enumerar todas as conexões únicas (gi, gj), o que sugere colocá-las todas em um contêiner associativo (não duplicado) e depois contar sua cardinalidade - em C ++, isso seria um std :: set <std :: pair <int, int>>. No seu pseudocódigo, você substituirá "nnz [i] ++" por "s.insert [pair (gi, gj)]" "e o número final de nonzeros será s.size (). Ele deve ser executado no tempo O (n-log-n), onde n é o número de não-zeros.

Como você provavelmente já conhece a variedade de possíveis soldados, você pode "exibir" a tabela pelo índice de soldados para melhorar o desempenho. Isso substitui o seu conjunto por um std :: vector <std :: set <int>>. Você preenche isso com "v [gi] .insert (gj)", então o número total de nonzeros vem da soma de v [gi] .size () para todos os gi's. Isso deve ser executado no tempo O (n-log-k), onde k é o número de incógnitas por elemento (seis para você - essencialmente uma constante para a maioria dos códigos pde, a menos que você esteja falando dos métodos hp).

(Nota - queria que este fosse um comentário sobre a resposta selecionada, mas demorou demais - desculpe!)

rchilton1980
fonte

0

Iniciar a partir da matriz esparsa $E^T$ de elementos de tamanho $\times$ DOFs.

E_{Eu j}^{T} = {\begin{array}{cc} 1 & Eu f d o f j \in e eu e m e n t Eu \\ 0 0 & e eu s e W h e r e \end{array}

$E_{ij}^T = \left\{ \begin{array}{cc} 1 & \mathrm{if\ dof}\ j \in \mathrm{element}\ i \\ 0 & \mathrm{elsewhere} \end{array} \right.$ Matriz

A = E E^{T}

$A = E\, E^T$ tem o padrão de esparsidade que você está procurando. Lembre-se de que implementar

E^{T}

$E^T$ é mais fácil, é por isso que eu defini

E^{T}

$E^T$ ao invés de

E

$E$ .

Nicola Cavallini
fonte

Cálculo da estrutura de escarsidade para matrizes de elementos finitos

Respostas: