Редактирование: MPI решение волнового уравнения

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

:<math>T'' = c\ddot{T}</math>

В конце программы добавлено однопроцессорное решение для проверки результатов вычислений.
<div class="mw-collapsible mw-collapsed" style="width:100%" >
'''Текст программы на языке C++ (разработчик [[Цветков Денис]]):''' <div class="mw-collapsible-content"> 
<syntaxhighlight lang="cpp" line start="1" enclose="div">
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <ctime>
#include "include/mpi.h"
using namespace std;

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

    // Объявление переменных
    int done = 0, n, myid, numprocs, i;
    int namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

    // Инициализация подсистемы MPI
    MPI_Init(&argc, &argv);
    // Получить размер коммуникатора MPI_COMM_WORLD
    // (общее число процессов в рамках задачи)
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
    // Получить номер текущего процесса в рамках
    // коммуникатора MPI_COMM_WORLD
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);
    MPI_Get_processor_name(processor_name,&namelen);

    int N = 20000 + 2;            // количество частиц в цепочке, 2 частицы для г.у.
    double t = 200;               // общее время расчета стержня
    double dt = 0.01;             // шаг

    // Измерение времени работы MPI
    double startwtime, endwtime, duration_MPI;
    if (myid == 0) {
        startwtime = MPI_Wtime();
    }

    // здесь задаются начальные условия
    double U[N], V[N], UU[N], VV[N];        // переменные UU[N], VV[N] для сбора конечных результатов
    for (int i = 1; i < N - 1; i++) {
        U[i] = 1;
        if (i < N / 4) V[i] = 0;
        else V[i] = 0.01;
    }

    int N_per_proc = ceil((N - 2) / numprocs);            // количество частиц на каждый процессор

    for (double tt = 0; tt < t; tt+= dt) {
        // зеркальные Г.У.
        U[0] = U[1];
        U[N-1] = U[N-2];

        // циклы с 1 частицы до N - 1, т.к. первая и последняя частицы используются для г.у.

        // расчет скоростей для данного шага
        for (int j = 1 + N_per_proc * myid; j <  1 + N_per_proc * (myid + 1); j++) {
            V[j] += (U[j + 1] - 2 * U[j] + U[j - 1]) * dt;
        }

        // расчет перемещений для данного шага
        for (int j = 1 + N_per_proc * myid; j <  1 + N_per_proc * (myid + 1); j++) {
            U[j] += V[j] * dt;
        }

        if (numprocs > 1)
        if (myid == 0) {
            MPI_Send(&U[N_per_proc], 1, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&U[N_per_proc + 1], 1, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        } else if (myid == (numprocs - 1)) {
            MPI_Recv(&U[N_per_proc * (numprocs - 1)], 1, MPI_DOUBLE, numprocs - 2, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            MPI_Send(&U[N_per_proc * (numprocs - 1) + 1], 1, MPI_DOUBLE, numprocs - 2, 0, MPI_COMM_WORLD);
        } else {
            MPI_Recv(&U[N_per_proc * myid], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            MPI_Send(&U[N_per_proc * myid + 1], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD);
            MPI_Send(&U[N_per_proc * (myid + 1)], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&U[N_per_proc * (myid + 1) + 1], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        }
	}

    MPI_Barrier(MPI_COMM_WORLD);            // здесь происходит синхронизация потоков, чтобы избежать состояния гонки
	MPI_Gather(&U[1 + N_per_proc * myid], N_per_proc, MPI_DOUBLE, &UU[1], N_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Gather(&V[1 + N_per_proc * myid], N_per_proc, MPI_DOUBLE, &VV[1], N_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Barrier(MPI_COMM_WORLD);            // здесь происходит синхронизация потоков, чтобы избежать состояния гонки

    // Измерение времени работы MPI
    if (myid == 0) {
        endwtime = MPI_Wtime();
        duration_MPI = (endwtime-startwtime);
    }

    // Если это главный процесс, вывод полученного результата
    if(myid==0) {

        // здесь происходит расчет системы на одном процессоре, для сравнения результатов
        clock_t start;
        double duration;
        start = clock();

        double U1[N], V1[N];
        for (int i = 1; i < N - 1; i++) {
            U1[i] = 1;
            if (i < N / 4) V1[i] = 0;
            else V1[i] = 0.01;
        }
        for (double i = 0; i < t; i+= dt) {
            U1[0] = U1[1];
            U1[N-1] = U1[N-2];

            for (int j = 1; j <  N - 1; j++) {
                V1[j] += (U1[j + 1] - 2 * U1[j] + U1[j - 1]) * dt;
            }

            for (int j = 1; j <  N - 1; j++) {
                U1[j] += V1[j] * dt;
            }
        }

        duration = ( clock() - start ) / (double) CLOCKS_PER_SEC;

        // вывод перемещений (однопроцессорное решение и MPI)
        for (int i = 1; i < N - 1; i++) {
            printf("i = %d, U1 = %f, U = %f\n", i, U1[i], UU[i]);
        }

        cout <<"duration: "<< duration <<"   duration_MPI: "<< duration_MPI <<'\n';
    }

    // Освобождение подсистемы MPI
    MPI_Finalize();
    return 0;
}
</syntaxhighlight>
</div>
</div>

== Сравнение времени работы программы ==
N = 20000;

t = 500;

dt = 0.01;

Ноутбук, двухядерный процессор (Intel Core i5-3317U CPU 1.70 GHz)

Один процессор: 8.409 с.
MPI: 4.986 c.
Ускорение ~69%
@@ Строка 138: / Строка 138: @@
 </div>
 </div>
-Результат работы программы (N = 100, t = 500, dt = 0.01):
-[[Файл: Rez mpi.png]]
-== Реализация распараллеливания ==
-Для того, чтобы каждый процесс мог производить вычисления параллельно другим процессам, стержень (состоящий из N частиц) делится на количество процессов. Каждый шаг вычисления процессы обмениваются данными о том, что происходит на границах вычисляемого ими участка стержня с помощью функций отправки (''MPI_Send'') и приема (''MPI_Recv'') данных.
-После окончания вычислений каждый процесс имеет у себя участок памяти с результатами вычислений принадлежащего процессу участка стержня. Для объединения этих данных используется функция ''MPI_Gather'', принцип работы данной функции отображен на рисунке ниже.
-[[Файл: Image003.jpg]]
-После сбора данных многопроцессорного вычисления они выводятся вместе с данными, вычисленными с помощью одного процесса, в сравнительной таблице.
-Чтобы сбор и вывод данных управляющим процессом (''root''-процессом) не начался до того, как остальные процессы завершат свои вычисления, используется функция барьерной синхронизации ''MPI_Barrier'', которая блокирует работу вызвавшего ее процесса до тех пор, пока все другие процессы группы также не вызовут эту функцию. Завершение работы этой функции возможно только всеми процессами одновременно.
-Для правильного расчета требуется, чтобы N - 2 было кратно количеству процессов (например, при 16 процессах можно использовать N = 32000 + 2).
 == Сравнение времени работы программы ==
+N = 20000;
-Время работы программы:
+t = 500;
-{| border=1 style="text-align:center"
-!bgcolor=#DDDDDD rowspan="2"|Процессор
-!bgcolor=#DDDDDD colspan="2"|Время, затраченное на частицу за шаг, нс
-!bgcolor=#DDDDDD rowspan="2"|~Ускорение, %
-|-
-!bgcolor=#DDDDDD|1 процесс
-!bgcolor=#DDDDDD|MPI
-|-
-!bgcolor=#DDDDDD|Intel Core i5-3317U CPU 1.70 GHz, 2 ядра
-|8.41
-|4.99
-!bgcolor=#FFFFBB|69
-|-
-!bgcolor=#DDDDDD|Суперкомпьютер, 48 ядер
-|3.02
-|0.131
-!bgcolor=#FFFFBB|2205
-|-
-|}
+dt = 0.01;
-График времени работы программы в зависимости от количества задействованных ядер суперкомпьютера:
+Ноутбук, двухядерный процессор (Intel Core i5-3317U CPU 1.70 GHz)
-[[Файл: T Nproc.png]]
+Один процессор: 8.409 с.
+MPI: 4.986 c.
+Ускорение ~69%