Решение двумерного уравнения Пуассона методом блочных итераций

• Введение
• 1. Задание
• 2. Теоретический материал
• 3. Реализация поставленной задачи
• 3.1 Блок-схема алгоритма
• 3.2 Параллельная программа
• 3.3 Анализ работы программы на разном числе процессоров
• 3.4 Коэффициент ускорения вычислений в зависимости от числа потоков
• 3.5 График изменения погрешности
• Выводы по работе
• Список использованной литературы

Решение двумерного уравнения Пуассона итерационным методом Зейделя

Необходимо найти численное решение задачи Дирихле для уравнения Пуассона

(1)

в прямоугольной области с граничными условиями

(2)

Для решения поставленной задачи нами была написана параллельная программа на языке программирования высокого уровня С с расширением MPI, а также аналогичная программа на OpenMP.

уравнение пуассон параллельная программа

1. Задание

Решение двумерного уравнения Пуассона итерационным методом Зейделя

Найти численное решение задачи Дирихле для уравнения Пуассона

(1)

в прямоугольной области с граничными условиями

(2)

1. Разработайте блок-схему реализации распараллеливания данного алгоритма и напишите параллельную программу на MPI для численного решения уравнения (1) с условиями (2) с помощью данной итерационной схемы. Используйте распараллеливание прогонки (встречная прогонка). Для хранения сеточной функции используйте два двумерных массива, целиком размещающихся в памяти процессоров. В одном из них размещайте (), во втором и после его обработки пересылайте все содержимое массива в предыдущий массив. Тем самым вы экономите память и имеете возможность считать до любого значения n.

2. На сетке 50×50 проведите расчеты на разном числе процессоров и постройте зависимость ускорения вычислений и затраты на межпроцессорные обмены в зависимости от числа процессоров. Найдите оптимальное соотношение между числом процессоров и ускорением счета. Эффективность параллельного алгоритма и его отладку следует проводить с использованием средств профилирования, разработанных на кафедре ВС СибГУТИ.

3. Напишите аналогичную программу на OpenMP, проведите расчета на сетке 50×50 и определите коэффициент ускорения вычислений в зависимости от числа потоков.

4. Постройте график изменения погрешности от числа итераций.

2. Теоретический материал

Блочный итерационный метод Зейделя

На равномерной прямоугольной сетке уравнение (1) аппроксимируется следующей разностной схемой

(4)

где, n — номер итерации.

Значения сеточной функции на границах области известно из граничных условий. Схему (4) можно записать в виде, удобном для реализации ее с помощью метода скалярной прогонки:

где

Значения прогоночных коэффициентов находятся по рекуррентным формулам, которые можно записать в виде:

.

Из граничных условий на левой границе определяются значения прогоночных коэффициентов .

После этого, учитывая, что, обратной прогонкой находятся все значения сеточной функции на n+1 — ом итерационном шаге:

Счет следует проводить прогонкой по оси ОХ (индекс i), начиная с индекса j = 1. В этом случае значение переменной известно из граничного условия. Окончанием итерационного процесса является выполнение условия

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

3. Реализация поставленной задачи

3.1 Блок-схема алгоритма

3.2 Параллельная программа

#include

#include

#include

#include

#define N1 50

#define N2 50

#define eps 0.1

double Y [N1 + 1] [N2 + 1], Ysh [N1 + 1] [N2 + 1];

double hx = 1. f / N1, hy = 2. f / N2;

/*Функия точного решения*/

double Fresh (double x, double y) {

return pow (x,

2) * pow (y,

2);

}

double RoFresh (double x, double y) {

return 2 * (pow (x,

2) + pow (y,

2));

}

/*Подпрограмма инициализации матрицы*/

void Inic () {

int i, j;

for (i = 0; i < N1 + 1; i++)

for (j = 0; j < N2 + 1; j++) {

if ((i! = 0) && (j! = 0) && (i! = N1) && (j! = N2))

Y [i] [j] = 0;

else

Y [i] [j] = Fresh ((i * hx), (j * hy));

}

}

int main (int argc, char **argv) {

int size, rank, flag = 1;

int i, j, f, it = 0;

double A = pow (hx,

2), B = A, D = pow (hy,

2), C = - 2. f / A — 2. f / D, F1, Fi, pogr = 0;

double t_c = 0.0, time = 0.0, s_t = 0.0;

double max, m;

double alfa [N-1], beta [N-1];

// double Y [N1 + 1] [N2 + 1], Ysh [N1 + 1] [N2 + 1];

// double hx = 1. f / N1, hy = 2. f / N2;

MPI_Status stat;

MPI_Init (&argc, &argv);

MPI_Comm_size (MPI_COMM_WORLD, &size);

MPI_Comm_rank (MPI_COMM_WORLD, &rank);

time — = MPI_Wtime ();

t_c — = MPI_Wtime ();

if (rank == 0) {

printf («%d n», size);

Inic ();

}

MPI_Barrier (MPI_COMM_WORLD);

MPI_Bcast (Y, (N1 + 1) * (N2 + 1), MPI_DOUBLE, 0, MPI_COMM_WORLD);

if (rank == 0) {

do {

for (i = 0; i <= N1; i++)

for (j = 0; j <= N2; j++)

Ysh [i] [j] = Y [i] [j];

for (i = 1; i <= (N1 — 1) / size; i++) {

for (j = 1; j < N2; j++) {

Fi = (- (Y [i+1] [j] + Y [i] [j+1])) / D + RoFresh;

}

}

for (i = 0; i <= (N — 1); i++) {

alfa += (- D [i]) / (C [i] + Ai*alfa [i-1]);

beta += (F [i] - A [i] *beta [i-1]) / (C [i] + Ai*alfa [i-1]);

}

for (j = (N+1); j <= 1; j—)

Ysh += alfa [i] *Y [i+1] [j] + beta [i];

if (size == 1) {

do {

for (i = 1; i < N1; i++) {

for (j = 1; j < N2; j++) {

Fi = (- (Y [i+1] [j] + Y [i] [j+1])) / D + RoFresh;

}

}

max = m = - 999;

for (i = 0; i <= N1; i++) {

for (j = 0; j <= N2; j++) {

pogr = fabs (Ysh [i] [j] - Y [i] [j]);

if (pogr > max) max = pogr;

pogr = fabs (Ysh [i] [j]);

if (pogr > m) m = pogr;

}

}

if (max / m < eps) {

t_c += MPI_Wtime ();

printf («Bce 4etKo!!! it = %d, time = %fn», it, t_c);

flag = 0;

}

// Перезапись данных

for (i = 1; i < N1; i++)

for (j = 1; j < N2; j++)

Y [i] [j] = Ysh [i] [j];

it++;

} while (flag);

exit (1);

} else {

s_t — = MPI_Wtime ();

// Посылка следующему процессору верхнюю строку из основной матрицы в нижнюю границу

MPI_Send (&Ysh, (N1 + 1) * (N2 + 1), MPI_DOUBLE, 1, 1, MPI_COMM_WORLD);

MPI_Recv (&Y, (N1 + 1) * (N2 + 1), MPI_DOUBLE, 1, 1, MPI_COMM_WORLD, &stat);

s_t += MPI_Wtime ();

}

} while (1);

}

if ((rank > 0) && (rank < size — 1)) {

do {

MPI_Recv (&Ysh, (N1 + 1) * (N2 + 1), MPI_DOUBLE, rank — 1, 1, MPI_COMM_WORLD, &stat);

for (i = (N1 — 1) / size * rank + 1; i < (N1 — 1) / size * (rank + 1) + 1; i++) {

for (j = 1; j < N2; j++) {

Fi = (- (Y [i+1] [j] + Y [i] [j+1])) / D + RoFresh;

}

}

s_t — = MPI_Wtime ();

MPI_Send (&Ysh, (N1 + 1) * (N2 + 1), MPI_DOUBLE, rank — 1, 1, MPI_COMM_WORLD);

MPI_Send (&Ysh, (N1 + 1) * (N2 + 1), MPI_DOUBLE, rank + 1, 1, MPI_COMM_WORLD);

MPI_Recv (&Y, (N1 + 1) * (N2 + 1), MPI_DOUBLE, rank + 1, 1, MPI_COMM_WORLD, &stat);

s_t += MPI_Wtime ();

} while (1);

}

if ((rank == size — 1) && (rank! = 0)) {

do {

s_t — = MPI_Wtime ();

MPI_Recv (&Ysh, (N1 + 1) * (N2 + 1), MPI_DOUBLE, rank — 1, 1, MPI_COMM_WORLD, &stat);

s_t += MPI_Wtime ();

// for (i = (N1−1) / size * (size — 1) + 1; i < N1; i++) {

for (i = (N1−1) / size * rank + 1; i < N1; i++) {

for (j = 1; j < N2; j++) {

Fi = (- (Y [i+1] [j] + Y [i] [j+1])) / D + RoFresh;

}

}

max = m = - 999;

for (i = 0; i <= N1; i++) {

for (j = 0; j <= N2; j++) {

pogr = fabs (Ysh [i] [j] - Y [i] [j]);

if (pogr > max) max = pogr;

pogr = fabs (Ysh [i] [j]);

if (pogr > m) m = pogr;

}

}

if (max / m < eps) {

time += MPI_Wtime ();

printf («Bce 4etKo!!! it = %d, time = %f, s_t = %f, time_s4eta = %f n», it, time, s_t, time-s_t);

exit (1);

}

s_t — = MPI_Wtime ();

MPI_Send (&Ysh, (N1 + 1) * (N2 + 1), MPI_DOUBLE, rank — 1, 1, MPI_COMM_WORLD);

s_t += MPI_Wtime ();

// Перезапись данных

for (i = 0; i <= N1; i++)

for (j = 0; j <= N2; j++)

Y [i] [j] = Ysh [i] [j];

it++;

} while (1);

}

MPI_Finalize ();

return 0;

}

3.3 Анализ работы программы на разном числе процессоров

Полученные результаты вычислений сведены в таблицу 1.

Таблица 1

На основе результатов вычислений был построен следующий график, изображенный на рисунке 1.

Рисунок 1 — Зависимость вычисления и передачи от числа процессоров

3.4 Коэффициент ускорения вычислений в зависимости от числа потоков

3.5 График изменения погрешности

Рисунок 2 — График изменения погрешности

Выводы по работе

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

При решении уравнения на 6 процессорах общее время вычисления составляет 0,63 107 условных единиц времени, из которых 6,43 527 тратится на пересылку данных между процессорами.

На основе результатов вычисления был построен график изменения погрешности от числа итерации. Из данного графика делаем вывод, что точность получаемый результатов зависит от кол-ва итераций — чем больше итераций, тем выше точность получаемых результатов.

1. Бахвалов Е. А., Жидков Н. П., Кобельков Г. Н. Численные методы: Учеб. пособие. — М: Наука, 1987. — 600 с.

2. Березин И. С., Жидков Н. П. Методы вычислений. т.1, т.2. — М.: Наука, 1997.

3. Волков Е. А. Численные методы: Учеб. пособие для вузов. М: Наука, 1987. — 248 с.

4. Воеводин В. В, Воеводин Вл.В. Параллельные вычисления. — Спб.: БХВ-Петербург, 2002. — 608 с.

5. Корнеев В. Д. Параллельное программирование в MPI. — Москва-Ижевск: Ин-т компьютерных исследований, 2003. — 304 с.