Параллельные вычисления с использованием нескольких ядер с Open-MP. - PullRequest
0 голосов
/ 09 сентября 2018

Я изо всех сил пытаюсь выяснить, как распараллелить этот код с OpenMP, любая помощь приветствуется. Ниже приведен базовый код и описание.

При моделировании совокупности мягких частиц (таких как белки в жидкости) существует сила отталкивания между парой частиц, когда они перекрываются. Цель этого задания - использовать параллельные вычисления для ускорения вычисления этих сил отталкивания, используя несколько ядер с Open-MP.

В функции силового отталкивания частицы имеют единичный радиус. Частицы находятся в «коробке моделирования» с размерами L × L × L. Размер L выбран так, что объемная доля частиц составляет φ = 0,3. Блок моделирования имеет периодические (циклические) граничные условия, что объясняет, почему нам нужно использовать функцию остатка для вычисления расстояния между двумя частицами. Если частицы перекрываются, то есть расстояние s между двумя частицами меньше 2, то сила отталкивания пропорциональна k (2-s), где k - постоянная силы. Сила находится вдоль вектора, соединяющего две частицы.

  1. Напишите программу, которая проверяет правильность вашего кода. Это может быть сделано путем вычисления правильных сил и сравнения их с силами, вычисленными вашим оптимизированным кодом. В своем отчете укажите, что ваша программа работает правильно, используя вашу тестовую программу
  2. Насколько быстрее ваш ускоренный код по сравнению с предоставленным базовым кодом? Включите время для разных размеров проблемы. Обязательно включите в свой отчет список своего кода.

Код для распараллеливания

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

double get_walltime() {
struct timeval tp;
gettimeofday(&tp, NULL);
return (double) (tp.tv_sec + tp.tv_usec*1e-6); }

void force_repulsion(int np, const double *pos, double L, double     krepulsion, double *forces)
{
int i, j;
 double posi [4]; double rvec [4];
 double s2, s, f;
 // initialize forces to zero
for (i=0; i<3*np; i++)
 forces[i] = 0.;
 // loop over all pairs
for (i=0; i<np; i++)
 {
posi[0] = pos[3*i ];
posi[1] = pos[3*i+1]; posi[2] = pos[3*i+2];
for (j=i+1; j<np; j++)
 {
// compute minimum image difference
 rvec[0] = remainder(posi[0] - pos[3*j ], L);
 rvec[1] = remainder(posi[1] - pos[3*j+1], L);
 rvec[2] = remainder(posi[2] - pos[3*j+2], L);
 s2 = rvec [0]* rvec [0] + rvec [1]* rvec [1] + rvec [2]* rvec [2];
 if (s2 < 4)
 {
 s = sqrt(s2);
 rvec[0] /= s; rvec[1] /= s;
 rvec[2] /= s;
 f = krepulsion*(2.-s);
 forces[3*i ] += f*rvec[0];
 forces[3*i+1] += f*rvec[1];
 forces[3*i+2] += f*rvec[2];
 forces[3*j ] += -f*rvec[0];
 forces[3*j+1] += -f*rvec[1];
 forces[3*j+2] += -f*rvec[2]; }
 } }
 }


 int main(int argc, char *argv[]) {
 int i;
 int np = 100; // default number of particles
 double phi = 0.3; // volume fraction
 double krepulsion = 125.; // force constant
 double *pos; double *forces;
 double L, time0 , time1;

if (argc > 1)
np = atoi(argv[1]);

L = pow(4./3.*3.1415926536*np/phi, 1./3.);
 // generate random particle positions inside simulation box
forces = (double *) malloc(3*np*sizeof(double));
pos = (double *) malloc(3*np*sizeof(double));
for (i=0; i<3*np; i++)
  pos[i] = rand()/(double)RAND_MAX*L;
// measure execution time of this function

time0 = get_walltime ();
force_repulsion(np, pos, L, krepulsion, forces);
time1 = get_walltime ();

printf("number of particles: %d\n", np);
printf("elapsed time: %f\n", time1-time0);
free(forces);
free(pos);
return 0; }

1 Ответ

0 голосов
/ 09 сентября 2018

Теоретически, это было бы так просто:

void force_repulsion(int np, const double *pos, double L, double krepulsion, 
 double *forces)
{

   // initialize forces to zero
   #pragma omp parallel for
   for (int i = 0; i < 3 * np; i++)
     forces[i] = 0.;
   // loop over all pairs
   #pragma omp parallel for
   for (int i = 0; i < np; i++)
   {
    double posi[4];
    double rvec[4];
    double s2, s, f;

    posi[0] = pos[3 * i];
    //...

Компиляция:

g++ -fopenmp example.cc -o example

Обратите внимание, что я не проверял правильность.Убедитесь, что у вас не будет глобальной переменной внутри параллели для (как я обновил ваш код ..)

...