Работа с разреженными матрицами в МКЭ

[Cyrille]

Дорогой народ!

Пишу свою конечно-элементную программу (на Фортране). Возник вопрос: когда я получаю локальные матрицы для каждого из элементов сетки, я их собираю в глобальные матрицы. Что лучше - сразу преобразовывать их к сжатому виду (имеется в виду компактное хранение), или сначала помещать в нормальные разреженные матрицы (предварительно разместив в памяти - allocate), потом преобразовывать их в сжатый вид, проводить нужные вычисления на данном шаге по времени уже в компактном виде, после чего удалять из памяти (deallocate) и снова-корова для нового шага по времени. Последнее, вообще говоря, выглядит ужасно. Но первый вариант пока не знаю, как реализовать. Кто с этим сталкивался, поделитесь опытом, пожалуйста!

Заранее благодарен всем отclickнувшимся.

[VLarin]

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

Когда добавляем локальные матрицы к глобальной - у нас есть локальные номера узлов и глобальные и мы можем заранее сказать, в какую позицию вектора, хранящего данные нашей sparse-матрицы, нужно добавить соответствующий элемент локальной матрицы.

Поэтому - 1. Нет необходимости в явном виде создавать все локальные матрицы и 2. Сборка может осуществляться прямо в матрицу в сжатом виде:)

[Cyrille]

Угум... Спасибо. А вот ещё вопрос. У меня схема - явная, и шаги по времени малы настолько, что свойства материалов фактически не меняются на большом интервале времени. Поэтому решать линейные уравнения на каждом шаге по времени невыгодно. Мне бы посчитать обратную матрицу (в моей задаче это матрица теплоёмкости), да просто домножать на неё правую часть на каждом шаге, обновляя её лишь изредка, но helas! - ни в одной библиотеке не вижу программ, обращающих разреженные матрицы в сжатом формате. Как быть, не подскажите? В принципе, можно сформировать обратную матрицу из решения СЛАУ с множественными правыми частями, но это опять будет несжатый формат ...

[VLarin]

Ну как, существует множество форм спарс-матриц, и, если я не ошибаюсь, есть формы представления, позволяющие легко вычислить обратную.
Здесь нужно просто литературу посмотреть - для начала можно Джоржа и Лью "Численное решение больших разреженных матриц".
Там и про формы представления матриц, и про перенумерацию узлов и прочее. Но она очень старая. Может быть есть и что-то более новое.

[Vassil]

В этом, как и в другие форумы, случается что обычно я (быть может и другие?) читаю только "свои" темы. А рядом в другой теме иногда может лежать нужный мне ответ ...

Для широкого круга задач, в том числе для компактного хранения редких матриц и векторов и для основных общих задач над ними, можете поискать процедуры например в Систематическом каталоге Библиотеки численного анализа НИВЦ МГУ:

http://www.srcc.msu.su/num_anal/lib_na/cat/cat0.htm

Посмотрите так же не помогут ли Вам

http://www.netlib.org/linalg/html_templ ... lates.html

?

[Comga]

Cyrille писал(а):

У меня схема - явная, и шаги по времени малы настолько, что свойства материалов фактически не меняются на большом интервале времени. Поэтому решать линейные уравнения на каждом шаге по времени невыгодно. Мне бы посчитать обратную матрицу (в моей задаче это матрица теплоёмкости), да просто домножать на неё правую часть на каждом шаге, обновляя её лишь изредка, но helas! - ни в одной библиотеке не вижу программ, обращающих разреженные матрицы в сжатом формате. Как быть, не подскажите? В принципе, можно сформировать обратную матрицу из решения СЛАУ с множественными правыми частями, но это опять будет несжатый формат ...

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

[Comga]

VLarin писал(а):

Здесь нужно просто литературу посмотреть - для начала можно Джоржа и Лью "Численное решение больших разреженных матриц".
Там и про формы представления матриц, и про перенумерацию узлов и прочее. Но она очень старая. Может быть есть и что-то более новое.

Да уж, очень старая. Посоветую Y.Saad. Iterative methods for sparse linear systems. Там же и пакет Sparskit или ITSOL.

[Cyrille]

Comga писал(а):

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

Итоговая система уравнений имеет вид
C*[T(j+1)-T(j)]/dtau = K*T(j)+F
Поэтому на каждом временном шаге приходится решать линейную систему C*T(j+1) = RHS.

Что обратная к разреженной будет полной, не подумал, спасибо - действительно, весь смысл теряется... То есть для больших матриц умножение на их обратные на отдельных шагах вместо решения СЛАУ с разреженными на каждом шаге не даст никакого эффекта?

Спасибо за наводки на книжки!

[Comga]

Cyrille писал(а):

Итоговая система уравнений имеет вид
C*[T(j+1)-T(j)]/dtau = K*T(j)+F
Поэтому на каждом временном шаге приходится решать линейную систему C*T(j+1) = RHS.

Что обратная к разреженной будет полной, не подумал, спасибо - действительно, весь смысл теряется... То есть для больших матриц умножение на их обратные на отдельных шагах вместо решения СЛАУ с разреженными на каждом шаге не даст никакого эффекта?

1. Что есть С? Матрица? Какой структуры? Не распадается ли она на блоки?
2. Уже для трехдиагональной матрицы обратная будет полной.
3. Эффект даст, отрицательный. Если будете грамотно решать систему.

[Cyrille]

C - это матрица NxN, где N - число узлов сетки. В данной задаче содержит значения теплоёмкости. Поскольку любой узел соединён не со всеми остальными, а только с несколькими узлами, матрица сильно разреженная. Разбивается ли она на блоки... мне кажется, нет, но про это я ещё почитаю, тут я плаваю немного.

Вообще, скорость счёта снижается именно из-за операций с большими матрицами, или из-за необходимости размещать в памяти такие громоздкие массивы? Или из-за неэффективного размещения полных разреженных матриц в памяти?

[Comga]

Cyrille писал(а):

Поскольку любой узел соединён не со всеми остальными, а только с несколькими узлами, матрица сильно разреженная

Все-таки непонятно, почему в явной схеме появляется эта матрица? Ведь в явной схеме все связи задействованы тольно на нижнем уровне!

[Cyrille]

Эта матрица появляется из дискретизации нестационарного уравнения теплопроводности
ro*Cp*dT/dtau = lambda*(d2T/dx^2 + d2T/dy^2) + qv
T(x) аппроксимируем с помощью метода конечных элементов, а Т(tau) - центральной конечной разностью. В итоге получаем систему C*dT/dtau = K*T + F.
В механике К называется матрицей жёсткости, в тепловых задачах это матрица теплопроводности, F - в механике вектор нагрузки, у нас - вектор (плотности) теплового потока. С - в механике забыл, как называется, у нас это матрица теплоёмкости.

[Comga]

Еще раз. В явных схемах нет связей на следующем уровне между значения в точках (только на текущем). Связь на следующем уровне есть только в неявных схемах. Собственно это отражено в их названиях.

[Аспир]

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

[Yuri Gendelman]

Аспир писал(а):

При этом матрица коэфициентов выходит плохо обусловленной.
Каким методом можна решать такую систему??

Решение можно сделать более устойчивым, если применять балансировку матрицы: http://www.srcc.msu.su/num_anal/lib_na/cat/cat529.htm.