Расчет разреженной структуры для матриц конечных элементов

13

Вопрос: Какие методы доступны для точного и эффективного расчета разреженной структуры матрицы конечных элементов?

Информация: я работаю над решателем уравнения давления Пуассона, использую метод Галеркина с квадратичной базой Лагранжа, написанный на C, и использую PETSc для хранения разреженных матриц и подпрограмм KSP. Чтобы эффективно использовать PETSc, мне нужно предварительно выделить память для глобальной матрицы жесткости.

В настоящее время я делаю пробную сборку, чтобы оценить количество ненулевых элементов в строке следующим образом (псевдокод)

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

Это, однако, переоценивает nnz, потому что некоторые взаимодействия узел-узел могут происходить в нескольких элементах.

Я подумал о том, чтобы попытаться отследить, какие взаимодействия i, j я обнаружил, но я не уверен, как это сделать без использования большого количества памяти. Я также мог бы зацикливаться на узлах и находить поддержку базисной функции в центре этого узла, но тогда мне пришлось бы искать все элементы для каждого узла, что кажется неэффективным.

Я нашел этот недавний вопрос, который содержал некоторую полезную информацию, особенно от Стефано М., который написал

Мой совет - реализовать его на python или C, применяя некоторые теоретические концепции графа, то есть рассматривать элементы в матрице как ребра графа и вычислять разреженную структуру матрицы смежности. Список списков или словарь ключей являются общими вариантами.

Я ищу больше деталей и ресурсов по этому вопросу. По общему признанию, я не знаю много теории графов, и я не знаком со всеми трюками CS, которые могли бы быть полезными (я подхожу к этому с математической стороны).

Благодарность!

matrix finite-element petsc performance Джон Эдвардсон
источник

5

Ваша идея отслеживать, какие i, j взаимодействия вы нашли, может сработать, я думаю, что это «трюк CS», о котором вы и Стефано М. имеете в виду. Это равносильно созданию вашей разреженной матрицы в формате списка списков .

Не знаю, сколько у вас CS, поэтому я извиняюсь, если это вам уже известно: в структуре данных связанного списка каждая запись хранит указатель на запись после нее и запись до. Дешево добавлять и удалять записи, но не так просто найти в них элементы - вам, возможно, придется просмотреть их все.

Итак, для каждого узла i вы храните связанный список. Затем вы перебираете все элементы; если вы обнаружите, что два узла i и j соединены, вы можете посмотреть в связанном списке i. Если j еще не существует, вы добавляете его в список, а также добавляете i в список j. Проще всего добавить их по порядку.

После того, как вы заполнили свой список списков, теперь вы знаете количество ненулевых записей в каждой строке матрицы: это длина списка этого узла. Эта информация - именно то, что вам нужно для предварительного выделения разреженной матрицы в структуре данных матрицы PETSc. Затем вы можете освободить свой список списков, потому что он вам больше не нужен.

Однако этот подход предполагает, что все, что у вас есть, это список узлов, которые содержит каждый элемент.

Некоторые пакеты генерации сетки - например, треугольник - могут выводить не только список элементов и какие узлы они содержат, но и список всех ребер в вашей триангуляции. В этом случае вы не рискуете переоценить количество ненулевых элементов: для кусочно-линейных элементов каждое ребро дает вам ровно 2 элемента матрицы жесткости. Вы используете кусочно-квадратичный, так что каждое ребро имеет значение для 4 записей, но вы поняли идею. В этом случае вы можете найти количество ненулевых записей в строке за один проход по краевому списку, используя обычный массив.

При таком подходе вы должны прочитать очень большой файл с жесткого диска, который на самом деле может быть медленнее, чем использование списка элементов, если ваши фактические вычисления не такие большие. Тем не менее, я думаю, что это проще.

Даниэль Шаперо
источник

Благодарю. У меня есть доступный список границ, поэтому я, вероятно, пока воспользуюсь вашим вторым методом, но я мог бы вернуться и попробовать первый метод, просто чтобы испачкать руки связанными списками и тому подобным (спасибо за вступление ... Я Мы взяли только базовый класс CS, и хотя я умею программировать, я не так много знаю о структурах данных и алгоритмах)

Джон Эдвардсон

Рад был помочь! Я получил много своих знаний по CS из этого: books.google.com/books?isbn=0262032937 - ради любви к Богу, читайте об амортизированном анализе. Программирование собственного связанного списка или структуры данных бинарного дерева поиска в C стоит того.

Даниэль Шаперо

5

Если вы укажете вашу сетку как DMPlex, а вашу разметку данных как PetscSection, то DMCreateMatrix () автоматически предоставит вам правильно предварительно выделенную матрицу. Вот примеры PETSc для проблемы Пуассона и Стокса .

Мэтт Кнепли
источник

2

Upvoted

Лично я не знаю ни одного дешевого способа сделать это, поэтому я просто переоцениваю число, т. Е. Использую достаточно большое значение для всех строк.

Например, для идеально структурированной сетки, состоящей из линейных шестнадцатеричных элементов с узлами, nnzs на строку как в диагональном, так и вне диагональном блоках имеют значение dof * 27. Для большинства полностью неструктурированных автоматически созданных шестнадцатеричных сеток это число редко превышает dof * 54. Для линейных тет у меня никогда не было необходимости выходить за пределы dof * 30. Для некоторых сеток с элементами очень плохой формы / с низким соотношением сторон вам, возможно, придется использовать немного большие значения.

Наказанием является то, что потребление локальной (по рангу) памяти составляет от 2x-5x, поэтому вам может потребоваться использовать больше вычислительных узлов в кластере, чем обычно.

Кстати, я попытался использовать списки с возможностью поиска, но время, необходимое для определения структуры разреженности, было больше, чем сборка / решение. Но моя реализация была очень простой и не использовала информацию о ребрах.

Другой вариант - использовать подпрограммы, такие как DMMeshCreateExodus, как показано в этом примере.

Stali
источник

0

Вы хотите перечислить все уникальные (gi, gj) соединения, которые предлагают поместить их все в (не дублирующий) ассоциативный контейнер и затем посчитать его мощность - в C ++ это будет std :: set <std :: pair <int, int>>. В вашем псевдокоде вы бы заменили «nnz [i] ++» на «s.insert [pair (gi, gj)]», и тогда окончательное число ненулевых элементов будет s.size (). Он должен работать за O (n-log-n) время, где n - число ненулевых.

Поскольку вы, вероятно, уже знаете диапазон возможных значений gi, вы можете «раскрутить» таблицу по индексу gi, чтобы повысить производительность. Это заменит ваш набор на std :: vector <std :: set <int>>. Вы заполняете это с помощью "v [gi] .insert (gj)", тогда общее число ненулевых значений получается из суммирования v [gi] .size () для всех gi. Это должно выполняться за время O (n-log-k), где k - количество неизвестных на элемент (шесть для вас - по сути, константа для большинства кодов pde, если вы не говорите о hp-методах).

(Примечание - хотел, чтобы это был комментарий к выбранному ответу, но это было слишком долго - извините!)

rchilton1980
источник

0

Начать с разреженной матрицы $E^T$ размерных элементов $\times$ DOFs.

Е_{я J}^{T} знак равно {\begin{array}{cc} 1 & я е d о е J \in е L е м е N T я \\ 0 & е L s е вес час е р е \end{array}

$E_{ij}^T = \left\{ \begin{array}{cc} 1 & \mathrm{if\ dof}\ j \in \mathrm{element}\ i \\ 0 & \mathrm{elsewhere} \end{array} \right.$ матрица

A = E E^{T}

$A = E\, E^T$ имеет шаблон разреженности, который вы ищете. Имейте в виду, что реализация

E^{T}

$E^T$ это проще, поэтому я определил

E^{T}

$E^T$ вместо того

E

$E$ ,

Никола Каваллини
источник

Расчет разреженной структуры для матриц конечных элементов

Ответы: