Адаптация алгоритма Верле под многопроцессорные системы

Компьютерное моделирование в современной физике решает ряд задач, которые непосильно разрешить постановкой натурного эксперимента. Современная физика описывает практически все закономерности макромира в нормальных условиях, а вести исследования в областях низких или высоких температур, давлений проблематично, т.к. нет аппаратуры, которая могла бы работать в таких условиях. Поэтому для решения ряда задач научное сообщество прибегает к численному эксперименту, который позволяет смоделировать с достаточной достоверностью физические процессы, протекающие в реальных системах и выявить необходимые закономерности [1].

Современные компьютеры обладают высокими вычислительными способностями, но для наиболее приближенного описания физических систем необходимо моделировать порядка 1023 частиц, что на данный момент времени выполнить невозможно. Использование алгоритмов параллельных вычислений позволяет существенно сократить время моделирования. Рассмотрим одну из возможностей адаптации алгоритма Верле моделирования методом молекулярной динамики (МД) под многопроцессорные вычислительные системы.

Метод молекулярной динамики, попросту говоря, является численной реализацией решения уравнений движения Ньютона для множества частиц. Полную силу, действующую на i -й атом, можно представить в форме суммы векторов [2]:

N -1

F i = - ∑ ∇Φ ( ij ) .                              (1)

При этом предполагается, что известны координаты центров всех ато-

мов и вид потенциала взаимодействия Ф ( r ). Энергия частиц инертного

газа, например, может определяться потенциалом Леннарда-Джонса:

Л" 112

Ф (r )=4^|_ r

<7

r

Если в некоторый момент времени t известны координаты и импульсы всех атомов, то с помощью уравнений Ньютона можно определить траекторию i -го атома на заданном промежутке времени. В случае отсутствия внешних полей его координата и скорость будут иметь вид

1 N - 1

qi ⁽ t + A t) = и ( t ) A t -- — ^ ^^{ф (}ij )( A t) ² + q i ⁽ t) ,             (3)

2 m_i j

1 N - 1

^u i ⁽ t + ^A t ) = ^u i ⁽ t) —Е^{УФ (} ij ^{) A} t t .                   ⁽⁴⁾

m i j

Возьмем для примера систему из N частиц инертного газа.

Алгоритм Верле позволяет повысить точность итерации моделирования за счет частичного изменения скоростей частиц [3]. В общем случае одна итерация данного алгоритма на языке C++ представляет вид: void Verlet(void){ for (int i=0; i!=N; i++){

/*Расчет координаты i-й частицы*/ particles[i].x = particles[i].x + particles[i].vx*dt + particles[i].ax*dt2_2;

particles[i].y = particles[i].y + particles[i].vy*dt + particles[i].ay* dt2_2;

particles[i].z = particles[i].z + particles[i].vz*dt + particles[i].az* dt2_2;

/*Корректировка координат частицы согласно условий периодических граничных условий*/ periodic(particles[i].x, particles[i].y, particles[i].z);

/*Частичное изменение скорости */ particles[i].vx+=particles[i].ax*dt*0.5;

particles[i].vy+=particles[i].ay*dt*0.5;

particles[i].vz+=particles[i].az*dt*0.5;

}

GetAccel();//расчет ускорений for (int i=0; i!=N; i++){

/*Изменение скорости с учетом вычисленного ускорения*/ particles[i].vx+=particles[i].ax*dt*0.5;

particles[i].vy+=particles[i].ay*dt*0.5;

particles[i].vz+=particles[i].az*dt*0.5;

}

}

В представленном листинге производится вызов подпрограммы GetAccel(), задачей которой является определение ускорений частиц ис- ходя из (4). Видно, что для определения сил, действующих на одну частицу, необходимо определять расстояния до всех остальных N-1 частиц. Эта процедура дает наибольшую вычислительную нагрузку. В классическом представлении на языке C++ код подпрограммы выглядит следующим образом:

void GetAccel(void){ for (int i=0; i!=N; i++){ particles[i].ax=0;

particles[i].ay=0;

particles[i].az=0;

} for (int i=0; i!=N-1; i++){ for (int j=i+1; j!=N; j++){ dx=particles[j].x-particles[i].x;

dy=particles[j].y-particles[i].y; dz=particles[j].z-particles[i].z; Separate(dx,dy,dz);

r=sqrt(dx*dx+dy*dy+dz*dz);

r=1/r; r2=r*r; r3=r2*r; r6=r3*r3;

a=24*r2*r6*(2*r6-1);

particles[i].ax+=a*dx;

particles[i].ay+=a*dy;

particles[i].az+=a*dz;

particles[j].ax-=a*dx;

particles[j].ay-=a*dy; particles[j].az-=a*dz;

}

}

}

В среде программирования MS Visual Studio имеется компонент Net. FrameworkBackgroundWorker, с помощью которого возможна организация вычислительных операций в асинхронном фоновом режиме.

Массив координат частиц системы при расчете ускорений представляет собой массив связанных переменных, здесь ускорение i -й частицы напрямую связано с положением любой другой j -й частицы, поэтому прямое разделение подпрограммы Get Accel () на два и более асинхронных процесса фактически невозможно. Вместе с тем возможно использование одного управляющего компонента BackgroundWorker0 и двух счетных компонентов BackgroundWorker1 и BackgroundWorker2. Тогда схема разделения вычислительных операций одной итерации МД будет представлять следующую последовательность:

• 1.    BackgroundWorker0 обнуляет ускорения на начале данной итерации, запускает BackgroundWorker1 и BackgroundWorker2 и ожидает событий завершения работы этих компонентов.

• 2.    С использованием компонента BackgroundWorker1 производится расчет ускорений по первым N /2 частицам системы, с использованием BackgroundWorker2 – по последним N /2 частицам.

• 3.    По завершении выполнения работы асинхронных потоков при помощи BackgroundWorker0 производится расчет ускорений, вызванных силами взаимодействия первых N /2 частиц и последних N /2 частиц.

Рис 1. Схема разделения подпрограммы GetAccel() на вычислительные потоки

Представленная схема использована для двухпроцессорной вычислительной системы. Теоретически прямой расчет ускорений для системы из тысячи частиц влечет за собой 499500 итераций цикла, при использовании предложенной схемы разделения вычислений на один процессор приходится 249500 итераций, что приводит к сокращению времени расчетов примерно в два раза. На практике увеличение скорости вычислений составляет порядка 60%. Это связано со спецификой работы операционной системы и неравнозначностью вычислительных ядер процессора.

Предложенная схема применима и для вычислительных комплексов с большим количеством процессоров. Для этого необходимо разделение вычислений на большее количество счетных потоков.