Опубликован: 11.02.2010 | Уровень: специалист | Доступ: платный | ВУЗ: Санкт-Петербургский государственный университет
Лекция 4:

Фортран 90

Лабораторная работа 3.1 Оптимизация и распараллеливание вычислительной программы на примере моделирования системы взаимодействующих частиц методом молекулярной динамики

Метод молекулярной динамики

Метод молекулярной динамики является одним из основных методов моделирования динамики систем, состоящих из взаимодействующих между собой частиц. Он применяется в молекулярной физике, астрофизике и других науках. Интерпретация понятия "частица" зависит от предметной области. В молекулярной физике это могут быть атом или молекула, а в астрофизике "частицами" являются звезды или планеты, взаимодействующие с другими астрономическими объектами посредством гравитационного взаимодействия. Динамика систем таких частиц часто описывается уравнениями классической механики.

Будем рассматривать систему N частиц, динамика которых описывается уравнением второго закона Ньютона:

F_i=m_i a_i

Здесь F_i - равнодействующая всех сил, действующих на i -ю частицу, m_i - ее масса и а_i - ускорение, с которым частица двигается под действием силы. Это обыкновенное дифференциальное уравнение второго порядка, которое можно записать в виде:

F_i=m_i\frac{d^2 r_i}{dt^2}

где r_i радиус-вектор i -й частицы.

Равнодействующая сил, является суммой парных взаимодействий i -й частицы со всеми остальными частицами:

F_i=F(r_i)=\sum\limits_{i\ne j=1}^N F(r_i,r_j)

или

F_i=\sum\limits_{j=1}^{i-1}F(r_i,r_j) + \sum\limits_{j=i+1}^N F(r_i,r_j)

Сила взаимодействия связана с потенциалом взаимодействия. В расчетах молекулярных систем часто используют потенциал Леннарда-Джонса:

V(r)=4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]

Здесь r =\mid r\mid , а \epsilon и \sigma - параметры потенциала. Первый определяет интенсивность взаимодействия, а второй - положение минимума (см. рис. 6.5).

Потенциал Леннарда-Джонса

Рис. 6.5. Потенциал Леннарда-Джонса

Потенциал Леннарда-Джонса хорошо описывает взаимодействие между атомами Аргона. Для аргона m\approx 6.69\times 10^{-23}\hat e \tilde a,\:\:\: \sigma\approx 3.405 \times 10^{-10} \grave i ,\:\:\: \epsilon \approx 1.654\times 10^{-21} \ddot A ae. Будем считать, что массы всех частиц одинаковы, система единиц выбрана таким образом, что m = 1, \epsilon=1 и \sigma=1. Тогда сила парного взаимодействия:

F(r_i,r_j)=\frac{24(r_i-r_j)}{{\mid r_i-r_j\mid }^2} \left[\frac{2}{{\mid r_i-r_j\mid }^{12}} - \frac{1}{{\mid r_i-r_j\mid }^6} \right]

Интерес могут представлять траектории частиц или термодинамические характеристики системы, например, зависимость температуры от времени:

T(t)=\frac{1}{3Nk_B}\sum\limits_{i=1}^N {\left \mid  \frac{dr_i(t)}{dt} \right \mid }^2

где k_B\approx 1.38\times 10^{-23}\frac{\ddot A ae}{\hat E} -постоянная Больцмана.

Для вычисления траекторий в методе молекулярной динамики используются различные алгоритмы. Они основаны на замене непрерывного времени дискретным набором значений t\to t^k=t^0+k\Delta t, \:\:\: k=1,2,\ldots.. Одним из алгоритмов является метод Верле:

r_i^{k+1}=2r_i^k-r_i^{k-1}+F(r_i^k)\frac{\Delta t^2}{m_i}

Здесь r_i^k=r_i(t^k). Наиболее трудоемким является вычисление силы, действующей на частицу.

Для оптимизации вычисления сил иногда вводится радиус обрезания - некоторый фиксированный (достаточно большой) радиус, такой, что пересчет сил, действующих между частицами, оказавшимися внутри этого радиуса, производится на каждом шаге. Пересчет сил, действующих со стороны более удаленных частиц, производится через несколько шагов. Такая стратегия позволяет уменьшить количество трудоемких операций расчета парных сил. Номера "ближайших" частиц могут храниться в специальном массиве. прТакой массив называется списком.Обычно используются два массива целого типа. Каждый элемент первого массива ( list_nb ) содержит количество ближайших соседей соответствующей частицы:

k0 = 0 
j0 = 0 
do i = 1,  n 
do j = 1, n
if(abs(x(i) - x(j)) < rcutoff)   then 
j0 = j0 + 1
list_nb(i)  = j end if 
end do
i_nb(i) = j0 - k0 end do

Во втором массиве содержатся номера соседей каждой частицы. Положения хранятся в элементах, начиная с start(i) и заканчивая finish(i):

start(1) = 1
finish(1) = i_nb(1) do i = 2, n
start(i) = finish(i - 1) + 1
finish(i)  = finish(i - 1) + i_nb(i) end do

Лабораторная работа

В заданиях лабораторной работы 3.1 предлагается выполнить оптимизацию вычислительной программы. Первый этап оптимизации - сокращение числа операций в последовательном коде на основе использования физических законов. Затем выполняется распараллеливание программы. Выбор инструмента распараллеливания должен быть аргументирован.

Данная работа не предполагает выполнения всех заданий. Учащемуся может быть предложен набор из нескольких заданий, обязательно включающий номера 1, 2 и 3, а также одно или несколько заданий, связанных с распараллеливанием.

По результатам составляется отчет, содержащий решения с обоснованием выбранного метода, результаты измерения быстродействия программы, а также физические результаты и выводы.

Выполнение данной работы рассчитано на несколько занятий (2-6 академических часов). Имеется полный исходный текст программы на языке Fortran 90.

Задания для практической работы
Задание 1

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

Задание 2

Откомпилировать программу, выполнить моделирование для системы с параметрами, заданными преподавателем. Определить процессорное время, потраченное на выполнение моделирования.

Задание 3

Выполнить оптимизацию последовательного кода программы на основе "физических соображений". Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 2.

Задание 4

Проанализировать последовательный код и выявить участки потенциального параллелизма. Выполнить распараллеливание с помощью OpenMP (если это возможно). Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 3.

Задание 5

Проанализировать последовательный код, выявить участки потенциального параллелизма и предложить способ распараллеливания в модели передачи сообщений. Выполнить распараллеливание с помощью MPI, используя блокирующие двухточечные операции обмена (если это возможно). Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатами, полученными в заданиях 3 и 4.

Задание 6

Проанализировать возможность использования других видов обмена (двухточечные неблокирующие, коллективные). Если предполагается выигрыш в производительности, модифицировать вариант программы с использованием MPI. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатами, полученными в заданиях 3 , 4 и 5.

Задание 7

Выполнить оптимизацию последовательного кода, полученного при выполнении задания 3, используя введение радиуса обрезания и метод генерации списков частиц. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 3.

Задание 8

Провести анализ зависимостей при генерации списков частиц. Выполнить распараллеливание, выбрав наиболее подходящий инструмент. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 7.

Задание 9

Выполнить распараллеливание вычисления сил, действующих на частицы, используя метод списков. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в заданиях 4 или 5.

Задание 10

Исследовать степень сбалансированности загрузки процессоров (ядер при использовании многоядерной архитектуры). Если сбалансированность неудовлетворительная, предложить способы ее улучшения.

Задание 11

Исследовать масштабируемость программы. Если масштабируемость неудовлетворительная, предложить способы ее улучшения.

Задание 12

На основании результатов, полученных при выполнении заданий данной лабораторной работы, написать отчет, в котором содержатся выводы об эффективности различных способов оптимизации исходного последовательного кода и трудоемкости реализации этих способов на практике.

Сергей Лебедев
Сергей Лебедев
Россия
Паулус Шеетекела
Паулус Шеетекела
Россия, ТГТУ, 2010