Obliczanie struktury rzadkości dla macierzy elementów skończonych

Pytanie: Jakie metody są dostępne w celu dokładnego i wydajnego obliczenia struktury rzadkości matrycy elementów skończonych?

Informacje: Pracuję nad solwerem Poissona Równania Ciśnienia, stosując metodę Galerkina z kwadratową podstawą Lagrange'a, napisaną w C, i używam PETSc do rzadkiego przechowywania macierzy i procedur KSP. Aby efektywnie korzystać z PETSc, muszę wstępnie przydzielić pamięć dla globalnej macierzy sztywności.

Obecnie wykonuję próbny zestaw, aby oszacować liczbę nonzeros na wiersz w następujący sposób (pseudokod)

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]++

Jednak przecenia to nnz, ponieważ niektóre interakcje węzeł-węzeł mogą wystąpić w wielu elementach.

Zastanawiałem się nad tym, aby sprawdzić, jakie interakcje i, j znalazłem, ale nie jestem pewien, jak to zrobić bez użycia dużej ilości pamięci. Mógłbym również zapętlić węzły i znaleźć wsparcie funkcji bazowej na środku tego węzła, ale wtedy musiałbym przeszukać wszystkie elementy dla każdego węzła, co wydaje się nieefektywne.

Znalazłem to ostatnie pytanie, które zawierało przydatne informacje, zwłaszcza od Stefano M., który napisał

radzę zaimplementować go w Pythonie lub C, stosując pewne teoretyczne koncepcje graficzne, tzn. rozważ elementy w macierzy jako krawędzie w grafie i oblicz strukturę rozproszeniową macierzy przylegania. Lista list lub słownik kluczy to często wybierane opcje.

Szukam więcej szczegółów i zasobów na ten temat. Wprawdzie nie znam zbyt wielu teorii grafów i nie znam wszystkich sztuczek CS, które mogą być przydatne (podchodzę do tego od strony matematycznej).

Dzięki!

Twój pomysł śledzenia, które interakcje i, j znalazłeś, może zadziałać, myślę, że to „sztuczka CS”, o której mówisz ty i Stefano M. Sprowadza się to do zbudowania rzadkiej macierzy w formacie listy list .

Nie jestem pewien, ile masz CS, więc przepraszam, jeśli jest to już znane: w strukturze danych z listą połączoną każdy wpis przechowuje wskaźnik do pozycji po niej i pozycji wcześniejszej. Dodawanie i usuwanie wpisów jest tanie, ale nie jest tak łatwe do znalezienia elementów - być może będziesz musiał przejrzeć je wszystkie.

Tak więc dla każdego węzła i przechowujesz połączoną listę. Następnie iterujesz przez wszystkie elementy; jeśli okaże się, że dwa węzły i i są połączone, przejrzyj listę połączoną i. Jeśli j jeszcze tam nie ma, dodajesz go do listy, a także dodajesz i do listy j. Najłatwiej jest je dodać po kolei.

Po zapełnieniu listy list znasz teraz liczbę niezerowych wpisów w każdym rzędzie macierzy: jest to długość listy tego węzła. Ta informacja jest dokładnie tym, czego potrzebujesz, aby wstępnie przydzielić rzadką macierz w strukturze danych macierzy PETSc. Następnie możesz zwolnić swoją listę list, ponieważ już jej nie potrzebujesz.

Jednak takie podejście zakłada, że ​​wszystko, co masz, to lista węzłów, które zawiera każdy element.

Niektóre pakiety generowania siatki - na przykład Trójkąt - mogą wyświetlać nie tylko listę elementów i zawartych w nich węzłów, ale także listę wszystkich krawędzi w twojej triangulacji. W takim przypadku nie ma ryzyka przeszacowania liczby niezerowych wpisów: w przypadku kawałków elementów liniowych każda krawędź daje dokładnie 2 wpisy macierzy sztywności. Używasz częściowego kwadratu, więc każda krawędź liczy się dla 4 wpisów, ale masz pomysł. W takim przypadku można znaleźć liczbę niezerowych wpisów na wiersz za pomocą jednego przejścia przez listę krawędzi przy użyciu zwykłej tablicy.

Przy takim podejściu musisz odczytać dodatkowy duży plik z dysku twardego, który może być wolniejszy niż użycie listy elementów, jeśli twoje obliczenia nie są tak duże. Niemniej jednak myślę, że to prostsze.

Jeśli określisz swoją siatkę jako DMPlex, a układ danych jako PetscSection, wówczas DMCreateMatrix () automatycznie dostarczy poprawnie wstępnie przydzieloną macierz. Oto przykłady PETSc dla problemu Poissona i problemu Stokesa .

Pozytywne

Osobiście nie znam żadnego taniego sposobu na zrobienie tego, więc po prostu przeceniam liczbę, tj. Używam rozsądnie dużej wartości dla wszystkich wierszy.

Np. Dla idealnie skonstruowanej siatki wykonanej z liniowych 8-węzłowych elementów sześciokątnych nnzs na rząd zarówno w blokach po przekątnej, jak i po przekątnej wynoszą dof * 27. Dla większości w pełni nieustrukturyzowanych automatycznie generowanych siatek heksadecymalnych liczba rzadko przekracza dof * 54. W przypadku tetów liniowych nigdy nie musiałem wychodzić poza dof * 30. W przypadku niektórych siatek z elementami o bardzo złym kształcie / niskim współczynniku kształtu może być konieczne użycie nieco większych wartości.

Kara polega na tym, że zużycie pamięci lokalnej (na poziomie) wynosi od 2x do 5x, więc może być konieczne użycie większej liczby węzłów obliczeniowych w klastrze niż zwykle.

Przy okazji próbowałem użyć przeszukiwalnych list, ale czas potrzebny na określenie struktury sparcytacji był większy niż montaż / rozwiązanie. Ale moja implementacja była bardzo prosta i nie wykorzystywała informacji o krawędziach.

Inną opcją jest użycie procedur takich jak DMMeshCreateExodus, jak pokazano w tym przykładzie.

Chcesz wyliczyć wszystkie unikalne połączenia (gi, gj), co sugeruje umieszczenie ich wszystkich w (nie powielającym się) pojemniku asocjacyjnym, a następnie policzenie jego liczności - w C ++ byłoby to std :: set <std :: pair <int, int>>. W swoim pseudokodzie zamieniłeś „nnz [i] ++” na „s.insert [pair (gi, gj)]”, a następnie końcową liczbą nonzeros jest s.size (). Powinien działać w czasie O (n-log-n), gdzie n jest liczbą nonzerów.

Ponieważ prawdopodobnie znasz już zakres możliwych Gi, możesz „rozstawić” tabelę za pomocą indeksu Gi, aby poprawić wydajność. To zastępuje twój zestaw std :: vector <std :: set <int>>. Wypełniasz to „v [gi] .insert (gj)”, wtedy całkowita liczba nonzeros pochodzi z sumowania v [gi] .size () dla wszystkich gi. Powinno to działać w czasie O (n-log-k), gdzie k jest liczbą niewiadomych przypadających na element (sześć dla ciebie - zasadniczo stała dla większości kodów pde, chyba że mówisz o metodach hp).

(Uwaga - chciałem, aby był to komentarz do wybranej odpowiedzi, ale był za długi - przepraszam!)

$E^T$$\times$