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

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

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

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

 1 #include <iostream>
 2 #include <stdio.h>
 3 #include <math.h>
 4 #include "include/mpi.h"
 5 using namespace std;
 6 
 7 int main(int argc, char *argv[]) {
 8 
 9     // Объявление переменных
10     int done = 0, n, myid, numprocs, i;
11     int namelen;
12     char processor_name[MPI_MAX_PROCESSOR_NAME];
13 
14     // Инициализация подсистемы MPI
15     MPI_Init(&argc, &argv);
16     // Получить размер коммуникатора MPI_COMM_WORLD
17     // (общее число процессов в рамках задачи)
18     MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
19     // Получить номер текущего процесса в рамках
20     // коммуникатора MPI_COMM_WORLD
21     MPI_Comm_rank(MPI_COMM_WORLD,&myid);
22     MPI_Get_processor_name(processor_name,&namelen);
23 
24     int N = 100 + 2;            // количество частиц в цепочке, 2 частицы для г.у.
25     double t = 30;              // общее время расчета стержня
26     double dt = 0.01;           // шаг
27 
28     // здесь задаются начальные условия
29     double U[N], V[N];
30     for (int i = 1; i < N - 1; i++) {
31         U[i] = 0;
32         if (i < N / 2) V[i] = 0;
33         else V[i] = 0.01;
34     }
35 
36     int N_per_proc = ceil(N / numprocs);            // количество частиц на каждый процессор
37 
38     // расчет системы
39     for (double i = 0; i < t; i+= dt) {
40         // зеркальные Г.У.
41         U[0] = U[1];
42         U[N-1] = U[N-2];
43 
44         // циклы с 1 частицы до N - 1, т.к. первая и последняя частицы используются для г.у.
45 
46         // расчет скоростей для данного шага
47         for (int j = 1 + N_per_proc * myid; j <  1 + N_per_proc * (myid + 1); j++) {
48             if (j < N - 1) {
49                 V[j] += (U[j + 1] - 2 * U[j] + U[j - 1]) * dt;
50             }
51         }
52 
53         MPI_Barrier(MPI_COMM_WORLD);            // здесь происходит синхронизация данных
54 
55         // расчет перемещений для данного шага
56         for (int j = 1 + N_per_proc * myid; j <  1 + N_per_proc * (myid + 1); j++) {
57             if (j < N - 1) {
58                 U[j] += V[j] * dt;
59             }
60         }
61 
62         MPI_Barrier(MPI_COMM_WORLD);            // здесь происходит синхронизация данных
63 	}
64 
65     // Если это главный процесс, вывод полученного результата
66     if(myid==0) {
67 
68         // здесь происходит расчет системы на одном процессоре, для сравнения результатов
69         double U1[N], V1[N];
70         for (int i = 1; i < N - 1; i++) {
71             U1[i] = 0;
72             if (i < N / 2) V1[i] = 0;
73             else V1[i] = 0.01;
74         }
75         for (double i = 0; i < t; i+= dt) {
76             U1[0] = U1[1];
77             U1[N-1] = U1[N-2];
78 
79             for (int j = 1; j <  N - 1; j++) {
80                 V1[j] += (U1[j + 1] - 2 * U1[j] + U1[j - 1]) * dt;
81             }
82 
83             for (int j = 1; j <  N - 1; j++) {
84                 U1[j] += V1[j] * dt;
85             }
86         }
87 
88         // вывод скорости узлов стержня (MPI и однопроцессорное решение)
89         for (int i = 1; i < N - 1; i++) {
90             printf("i = %d, V = %f, V1 = %f\n", i, V[i], V1[i]);
91         }
92     }
93 
94     // Освобождение подсистемы MPI
95     MPI_Finalize();
96     return 0;
97 }