MPI решение волнового уравнения
Здесь представлено стохастическое решение волнового уравнения:
В конце программы добавлено однопроцессорное решение для проверки результатов вычислений.
1 #include <iostream>
2 #include <stdio.h>
3 #include <math.h>
4 #include <ctime>
5 #include "include/mpi.h"
6 using namespace std;
7
8 int main(int argc, char *argv[]) {
9
10 // Объявление переменных
11 int done = 0, n, myid, numprocs, i;
12 int namelen;
13 char processor_name[MPI_MAX_PROCESSOR_NAME];
14
15 // Инициализация подсистемы MPI
16 MPI_Init(&argc, &argv);
17 // Получить размер коммуникатора MPI_COMM_WORLD
18 // (общее число процессов в рамках задачи)
19 MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
20 // Получить номер текущего процесса в рамках
21 // коммуникатора MPI_COMM_WORLD
22 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
23 MPI_Get_processor_name(processor_name,&namelen);
24
25 int N = 20000 + 2; // количество частиц в цепочке, 2 частицы для г.у.
26 double t = 200; // общее время расчета стержня
27 double dt = 0.01; // шаг
28
29 // Измерение времени работы MPI
30 double startwtime, endwtime, duration_MPI;
31 if (myid == 0) {
32 startwtime = MPI_Wtime();
33 }
34
35 // здесь задаются начальные условия
36 double U[N], V[N], UU[N], VV[N]; // переменные UU[N], VV[N] для сбора конечных результатов
37 for (int i = 1; i < N - 1; i++) {
38 U[i] = 1;
39 if (i < N / 4) V[i] = 0;
40 else V[i] = 0.01;
41 }
42
43 int N_per_proc = ceil((N - 2) / numprocs); // количество частиц на каждый процессор
44
45 for (double tt = 0; tt < t; tt+= dt) {
46 // зеркальные Г.У.
47 U[0] = U[1];
48 U[N-1] = U[N-2];
49
50 // циклы с 1 частицы до N - 1, т.к. первая и последняя частицы используются для г.у.
51
52 // расчет скоростей для данного шага
53 for (int j = 1 + N_per_proc * myid; j < 1 + N_per_proc * (myid + 1); j++) {
54 V[j] += (U[j + 1] - 2 * U[j] + U[j - 1]) * dt;
55 }
56
57 // расчет перемещений для данного шага
58 for (int j = 1 + N_per_proc * myid; j < 1 + N_per_proc * (myid + 1); j++) {
59 U[j] += V[j] * dt;
60 }
61
62 if (numprocs > 1)
63 if (myid == 0) {
64 MPI_Send(&U[N_per_proc], 1, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
65 MPI_Recv(&U[N_per_proc + 1], 1, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
66 } else if (myid == (numprocs - 1)) {
67 MPI_Recv(&U[N_per_proc * (numprocs - 1)], 1, MPI_DOUBLE, numprocs - 2, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
68 MPI_Send(&U[N_per_proc * (numprocs - 1) + 1], 1, MPI_DOUBLE, numprocs - 2, 0, MPI_COMM_WORLD);
69 } else {
70 MPI_Recv(&U[N_per_proc * myid], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
71 MPI_Send(&U[N_per_proc * myid + 1], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD);
72 MPI_Send(&U[N_per_proc * (myid + 1)], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD);
73 MPI_Recv(&U[N_per_proc * (myid + 1) + 1], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
74 }
75 }
76
77 MPI_Barrier(MPI_COMM_WORLD); // здесь происходит синхронизация потоков, чтобы избежать состояния гонки
78 MPI_Gather(&U[1 + N_per_proc * myid], N_per_proc, MPI_DOUBLE, &UU[1], N_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD);
79 MPI_Gather(&V[1 + N_per_proc * myid], N_per_proc, MPI_DOUBLE, &VV[1], N_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD);
80 MPI_Barrier(MPI_COMM_WORLD); // здесь происходит синхронизация потоков, чтобы избежать состояния гонки
81
82 // Измерение времени работы MPI
83 if (myid == 0) {
84 endwtime = MPI_Wtime();
85 duration_MPI = (endwtime-startwtime);
86 }
87
88 // Если это главный процесс, вывод полученного результата
89 if(myid==0) {
90
91 // здесь происходит расчет системы на одном процессоре, для сравнения результатов
92 clock_t start;
93 double duration;
94 start = clock();
95
96 double U1[N], V1[N];
97 for (int i = 1; i < N - 1; i++) {
98 U1[i] = 1;
99 if (i < N / 4) V1[i] = 0;
100 else V1[i] = 0.01;
101 }
102 for (double i = 0; i < t; i+= dt) {
103 U1[0] = U1[1];
104 U1[N-1] = U1[N-2];
105
106 for (int j = 1; j < N - 1; j++) {
107 V1[j] += (U1[j + 1] - 2 * U1[j] + U1[j - 1]) * dt;
108 }
109
110 for (int j = 1; j < N - 1; j++) {
111 U1[j] += V1[j] * dt;
112 }
113 }
114
115 duration = ( clock() - start ) / (double) CLOCKS_PER_SEC;
116
117 // вывод перемещений (однопроцессорное решение и MPI)
118 for (int i = 1; i < N - 1; i++) {
119 printf("i = %d, U1 = %f, U = %f\n", i, U1[i], UU[i]);
120 }
121
122 cout <<"duration: "<< duration <<" duration_MPI: "<< duration_MPI <<'\n';
123 }
124
125 // Освобождение подсистемы MPI
126 MPI_Finalize();
127 return 0;
128 }
Результат работы программы (N = 100, t = 500, dt = 0.01):
Реализация распараллеливания[править]
Для того, чтобы каждый процесс мог производить вычисления параллельно другим процессам, стержень (состоящий из N частиц) делится на количество процессов. Каждый шаг вычисления процессы обмениваются данными о том, что происходит на границах вычисляемого ими участка стержня с помощью функций отправки (MPI_Send) и приема (MPI_Recv) данных.
После окончания вычислений каждый процесс имеет у себя участок памяти с результатами вычислений принадлежащего процессу участка стержня. Для объединения этих данных используется функция MPI_Gather, принцип работы данной функции отображен на рисунке ниже.
После сбора данных многопроцессорного вычисления они выводятся вместе с данными, вычисленными с помощью одного процесса, в сравнительной таблице.
Чтобы сбор и вывод данных управляющим процессом (root-процессом) не начался до того, как остальные процессы завершат свои вычисления, используется функция барьерной синхронизации MPI_Barrier, которая блокирует работу вызвавшего ее процесса до тех пор, пока все другие процессы группы также не вызовут эту функцию. Завершение работы этой функции возможно только всеми процессами одновременно.
Для правильного расчета требуется, чтобы N - 2 было кратно количеству процессов (например, при 16 процессах можно использовать N = 32000 + 2).
Сравнение времени работы программы[править]
Время работы программы:
| Процессор | Время, затраченное на частицу за шаг, нс | ~Ускорение, % | |
|---|---|---|---|
| 1 процесс | MPI | ||
| Intel Core i5-3317U CPU 1.70 GHz, 2 ядра | 8.41 | 4.99 | 69 |
| Суперкомпьютер, 48 ядер | 3.02 | 0.131 | 2205 |
График времени работы программы в зависимости от количества задействованных ядер суперкомпьютера:
