MPI решение волнового уравнения

Здесь представлено стохастическое решение волнового уравнения:

[math]T'' = c\ddot{T}[/math]

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

Текст программы на языке C++ (разработчик Цветков Денис):

  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	t	dt
20000	500	0.01

Время работы программы:

Процессор	1 процесс	MPI	~Ускорение, %
Intel Core i5-3317U CPU 1.70 GHz, 2 ядра	8.409 с	4.986 с	69